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Abstract 

The  goal  of  this  project  was  to  develop  a  set  of  analysis  tools  for  performance  prediction  of 
flexural-disk  acoustic  projectors.  The  emphasis  is  on  analytically  based  formulations  designed  to 
merge  the  physics  of  the  flexural-disk  structure  with  typical  operational  requirements  of 
frequency,  source  level,  bandwidth,  and  operating  depth.  While  finite-element  analysis  coupled 
with  acoustic  radiation  models  would  be  used  to  refine  designs,  analytical  models  permit 
isolation  of  the  major  effects  of  essential  design  variables.  By  developing  an  understanding  of 
the  transducer  performance,  the  effectiveness  of  various  strategies  for  flexural-disk  designs  can 
be  evaluated. 

The  model  described  here  covers  four  configurations  of  flexural-disk  transducer  element.  These 
are  two-  and  three-layer  disks  clamped  at  the  center  and  two-  and  three-layer  disks  simply 
supported  at  the  edge.  One  of  the  layers  in  each  case  is  a  passive  substrate,  the  other  layer(s)  is 
(are)  piezoelectric  ceramic.  The  only  constraints  on  the  layer  thickness  are  that  the  three-layer 
configurations  must  have  the  same  thickness  of  ceramic  on  both  sides  of  the  substrate  and  that  in 
all  configurations  the  overall  thickness  be  small  with  respect  to  the  radius.  For  the  center- 
clamped  case,  a  non-zero  radius  for  the  clamping  post  is  permitted  (required,  in  fact,  since  the 
zero-radius  center  clamp  is  unrealistic).  An  option  is  included  to  permit  treatment  of  a  fluid- 
filled  volume  behind  the  disk. 

The  equivalent-circuit  model  is  constructed  including  dielectric  and  mechanical  losses  and  it  also 
includes  a  simple  model  for  the  radiation  load.  Using  this  equivalent-circuit  model,  the  driving- 
point  admittance  (with  and  without  water  load),  the  transmitting  voltage  and  current  responses, 
and  the  free-field  voltage  receiving  response  are  calculated.  With  the  file  structure  employed  in 
the  MathCad  program,  the  user  can  export  the  equivalent-circuit  parameters  to  other  applications 
or  the  user  can  import  equivalent-circuit  parameters  from  other  sources  to  run  the  internal 
performance  calculations. 
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I.  Introduction 

This  report  describes  the  development  of  an  equivalent  circuit  model  for  flexural  disk 
transducers'.  Four  configurations  are  considered  -  two  center-supported  structures  (Figs.  1  and 
2)  and  two  edge-supported  structures  (Figs.  3  and  4). 


Fig.  1 .  Two-layer,  center-supported  flexural  disk  structure.  In  this  view,  the  water  would  be  on 
the  upper  face  and  the  lower  face  would  contact  a  gas  back-volume.  For  many  practical 
transducers,  the  structure  would  be  doubled  so  as  to  have  two  structures  back-to-back.  The  edges 
are  nominally  free  but  would,  in  practice,  be  sealed  with  boots  or  bellows. 


Fig.  2.  Three-layer,  center-supported  flexural  disk  structure.  This  structure  is  assumed  to  be 
symmetric  with  equal  thickness  of  ceramic  on  either  side  of  the  passive  substrate. 
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SUPPORT 

Fig.  3.  Two-layer,  edge-supported  flexural  disk  structure.  As  in  Fig.  1,  water  would  contact  the 
upper  surface.  For  designs  not  compensated  for  hydrostatic  pressure,  this  would  ensure  that  the 
ceramic  is  in  compression.  For  the  edge-supported  transducers,  the  housing  contacts  the  active 
element  at  its  circumference.  The  edge  is  simply  supported. 


HINGED  EDGE 
SUPPORT 


Fig.  4.  Three-layer,  edge-supported  flexural  disk  structure.  As  in  Fig.  2,  only  the  symmetric  case 
with  both  ceramic  disks  equal  in  thickness  will  be  considered. 

Matrix  notation  is  used  throughout.  In  order  to  provide  continuity  with  the  bulk  of  the 
transducers  literature,  electric  flux  density,  D,  is  used  to  combine  the  effects  of  free  space  and 
materials  properties  (through  the  permittivity,  e).  A  better  approach  is  to  treat  all  dielectric 
properties  through  the  polarization,  P,  and  susceptibility,  %■  As  long  as  the  material  is  strictly 
linear,  use  of  the  older  convention  will  not  cause  problems  but,  for  high  drive  levels,  a 
formulation  in  terms  of  polarization  may  be  of  more  value2,3.  Also,  thin-plate  theory4  will  be 
used;  however,  the  approach  to  including  shear  will  be  described  briefly. 


II.  Basic  Piezoelectric  Equations 

To  solve  the  flexural-disk  problem,  the  plane-stress  assumptions  are  invoked  and  the  matrix 
equations  can  be  reduced  significantly.  The  basic  assumptions  of  plane-stress  are: 

•  Normal  stress  in  poling  direction  ( Tz  or  F)  is  negligible. 

•  “Vertical”  shear  stresses  {Ta,  Ty*  or  f4,  Ts)  are  negligible. 
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Furthermore,  from  the  electrode  configuration  (electrodes  only  on  disk  surfaces,  not  edges),  the 
electric  field  is  assumed  to  be  significant  only  in  poling  direction  (Ez  or  £3). 


The  basic  equations  that  couple  the  mechanical  and  piezoelectric  behavior  relate  four  variables: 
the  strain,  S,  the  stress,  T,  the  electric  field,  E,  and  the  electric  flux  density,  D.  In  three 
dimensions,  each  of  these  quantities  is  a  matrix.  Since  the  plane-stress  assumptions  constrain 
stress  not  strain,  the  following  general  form5  is  the  appropriate  starting  point: 


S  =  sE  T  +  d,r  E 


(1) 


D  =  dT  +  eT  E 


(2) 


which  can  be  expanded  without  further  approximation  as 


Si 

=  sEJ, 

+ 

4t2 

+ 

dlx  E3 

(3) 

s2 

=  4^ 

+ 

4j2 

+ 

Aei 

(4) 

A 

=  dlxTx 

+ 

4J2 

+ 

eT  E 

fc3j  U 3 

(5) 

This  is  a  complete  description  of  the  in-plane  mechanical  behavior  and  the  electrical  behavior. 
The  out-of-plane  strain,  £3,  is  not  zero  but  that  strain  does  not  enter  our  analysis.  If  we  had 
chosen  an  equation  set  with  5  as  an  independent  variable,  then  S3  could  be  eliminated  using  the 
£3  =  0  equation. 

In  some  circumstances,  other  forms  of  these  equations  are  more  convenient.  If  we  treat  Eqs.  3-5 
as  a  single  matrix  equation6, 


'sx' 

'4 

4 

< 

'A 

S2 

— 

4_ 

4 

4. 

A. 

k 

XI 

4j 

A. 

(6) 


where  the  3x3  matrix  on  the  right-hand  side  has  been  partitioned  to  show  its  composition. 
Inverting  this  matrix  equation  leads  directly  to  the  equivalent  coefficients  for  the  plane-stress 
approximation: 


'A 

c4 

ccD 

cc12 

-hK 

'5,' 

t2 

— 

Az 

cc,i 

-Mh> 

A. 

rhh 31 

-hK 

Tub 

1 _ 

A. 

The  effective  stiffness  coefficients  ( cc ),  piezoelectric  coefficient  (hh),  and  inverse  permittivity 
(PP)  can  be  extracted  directly  from  the  inverse  of  the  original  3x3  matrix.  In  order  to  emphasize 
the  difference  between  the  two-  and  three-dimensional  coefficients,  those  coefficients  having 
different  values  in  two  dimensions  will  be  denoted  by  doubling  the  letter.  For  example,  the 
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stiffness  coefficients  for  the  reduced  two-dimensional  problem  will  be  denoted  cci  i,  and  cc\z.  In 
setting  up  these  disk  problems,  manufacturers’  values  for  the  coefficients  can  only  be  used  in  Eq. 
3.  The  other  coefficients  must  be  derived  through  the  matrix  inverse.  The  matrix  equation,  Eq. 
7,  contains  the  three  equations, 


Tx  =  ccxDx  5,  +  cc°  S2  -  hh^  D2 

T2  =  ccX2  5,  +  ccfl  S2  -  hh3l  D} 

E,  =  -A/t,,  S',  -  hh^S,  +  pfcD, 


Also,  in  polar  coordinates,  the  strain  matrix  is 


with  an  equivalent  form  for  the  stress  matrix. 
Summary  of  Conversion  to  Two-Dimensions 


(8) 

(9) 

(10) 


(ID 


In  order  to  construct  the  plane-stress  equations,  the  following  subset  of  the  three-dimensional 
quantities  are  required: 

•s£ii>  s£i2»  dzu  and  e7^. 


These  permit  construction  of  the  matrix  equation,  Eq.  6.  Inverting  the  3x3  matrix  in  Eq.  6 
produces  the  values  for  ccD\\,  ccD\i,  hh^,  and  Pffa. 

To  illustrate  the  difference  between  the  effective  plane-stress  quantities  and  the  three- 
dimensional  quantities,  consider  PZT-8.  Table  I  compares  the  2-D  and  3-D  values  for  the 
stiffnesses,  the  hu  coefficient,  and  the  permittivity  at  constant  strain. 

Table  I.  Comparison  of  properties7  for  three-dimensional  analysis  and  two- 
dimensional  (plane-stress)  analysis  for  PZT8. 


ProDertv 

3-D 

2-D 

lunitsl 

11.5 

same 

X  10  A2m2/N 

/,2 

-3.7 

same 

x  10 A2m2/N 

<*3. 

-97 

same 

x  WnC/N 

£ 33/69 

1000 

same 

cD  \\ 

147 

97 

x  10 9  N/m2 

c\  2 

86.8 

55.2 

x  IQ9  N/m2 

A31 

-0.78 

-1.93 

x  109  V/m 

6*33/69 

561 

727 

- 
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A  point  to  note  regarding  two-dimensional  analysis  of  piezoelectric  materials  concerns  the 
Poisson’s  ratio,  a.  There  are  two  values  for  a  -  the  value  for  constant  electric  field,  (P,  and  the 
value  for  constant  electric  flux  density  (or  constant  polarization),  cP .  For  either  value,  the 
Poisson’s  ratio  in  the  12-plane  is 


(12) 


Because  the  piezoelectric  material  is  anisotropic,  the  Poisson’s  ratio  in  the  12-plane  can  be  larger 
than  0.5.  The  Poisson’s  ratio  in  either  the  13-  or  23-plane  is  considerably  smaller  so  the  material 
is  still  volumetrically  stable.  For  PZT-8  (see  Table  I),  the  value  of  cP  derived  from  the  3-D 
properties  is  about  0.55.  The  corresponding  value  from  the  2-D  properties  is  0.46  but,  for  some 
materials,  the  cP  from  the  2-D  properties  is  also  greater  than  0.5. 


Table  II  presents  a  summary  of  the  appropriate  two-dimensional  properties  for  several  generic 
piezoceramic  types7. 


Table  II.  Properties  for  typical  piezoceramic  materials.  Properties  with  doubled  letters  are  those 
values  corresponding  to  plane-stress  conditions. 


Property 

PZT4 

PZT5H 

PZT8 

lunitsl 

12.3 

16.5 

11.5 

x  10  'i2m2/N 

/.2 

-4.05 

-4.78 

-3.7 

x  10‘12  m2/N 

d. „ 

-122 

-274 

-97 

x  10  'nm2/N 

Pi-i/eo 

1300 

3400 

1000 

- 

P 

7600 

7500 

7500 

kg/m 3 

SD  U 

11.0 

14.0 

10.4 

x  10’12  m2/N 

P\2 

-5.34 

-7.27 

-4.76 

x  10'12  m2/N 

CC°  |i 

119 

97.8 

121 

x  10 9  N/m2 

CCD\1 

57.7 

50.8 

55.2 

x  10 9  N/m2 

hhji 

-1.87 

-1.35 

-1.93 

x  109  V/m 

eP^/Eq 

892 

1950 

728 

- 

III.  Relationship  of  Strain  to  Deflection 

In  order  to  analyze  the  behavior  of  flexural-disk  structures,  a  deflection  function  is  determined. 
For  thin-plate  problems,  the  exact  solution  can  be  found  for  both  static  deflection  and  for 
dynamic  deflection  at  resonance.  Alternatively,  approximate  forms  can  be  developed  for  the 
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deflection  function  and  optimized  by  minimizing  the  strain  energy  (for  static  deflection)  or  by 
minimizing  the  resonance  frequency  (for  resonance  deflection).  Once  the  deflection  function  is 
determined,  the  strains  can  be  evaluated. 

In  polar  coordinates,  if  the  deflection  function  is  w(r,  9),  then 


il 

= 

(d2w) 

■*W 

(13) 

s2  = 

^96  ~ 

—z 

Idw  1 

rdr  +  r2 

d2w ) 

W’J 

(14) 

*6  = 

& 

II 

-2  z 

(\  a2w 
rdrae 

1 

T2'dG^ 

(15) 

where  z  is  the  displacement  (in  the  3-direction)  away  from  the  neutral  plane. 

In  general,  flexural-disk  transducers  will  exploit  only  axially  symmetric  modes.  In  those  cases, 
the  strains  are: 


S2 


Jee 


—z 


a2  v»\ 

dr2 

V  J 


1  3w 
r  dr 

J 


(16) 


(17) 


and  S6  is  zero. 
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IV.  Basic  Solution  Process 


The  analysis  starts  with  determination  of  the  deflection  function.  For  simple,  homogeneous  disk 
problems  an  exact  solution8  for  deflection  can  be  done  (see  Appendix  A).  The  Version  3  model 
uses  an  approximate  technique  (see  Appendix  B)  to  permit  flexibility  in  future  revisions.  Once 
the  deflection  function  has  been  determined,  the  strain  is  derived  from  Eqs.  16  and  17.  Then  the 
basic  piezoelectric  equation  set  with  strain  as  an  independent  variable  (Eqs.  8-10)  can  be  solved. 


Solution  of  the  basic  piezoelectric  equation  set  proceeds  as  follows: 


A.  Integrate  (10)  over  z  (thickness)  and  solve  for  £>3.  Because  there  are  no  electrodes  in  the 
z-direction  (i.e.,  on  the  edges  of  the  plate),  £>3  is  not  a  function  of  z.  The  electric  field, IE, 
is  a  function  of  z;  however,  it  is  not  necessary  to  determine  that  dependence.  The  integral 
of  E  over  thickness  is  identical  to  the  electrode  voltage,  which  is  the  quantity  of 
importance. 

B.  Integrate  this  expression  for  £>3  over  the  electrode  area.  The  electrode  voltage  is  constant 
over  the  electrode.  The  integral  of  D  over  the  electrode  area  is  the  total  charge.  Set  the 
charge  equal  to  zero  and  solve  for  the  open-circuit  voltage,  V°c. 

C.  The  result  produced  in  A  gives  D3  as  a  function  of  electrode  voltage  and  the  strains.  D^oc 
results  from  setting  V  =  V°c\  D/c  results  from  setting  V=0. 

D.  The  potential  energy  density  for  strain  energy  is  Vi  S\  T\  +  Vi  S2T2.  The  open-circuit 
potential  energy  density  is  written  using  (8)  and  (9)  with  D°c\  the  short-circuit  potential 
energy  density  is  written  with  DiC . 

E.  The  strains,  Si  and  S2,  can  be  written  in  terms  of  the  deflection  function  according  to  Eqs. 
16  and  17. 

F.  Integrate  the  potential  energy  densities  over  the  ceramic  volume  to  obtain  U°c  and  Usc. 

G.  Integrate  D}SC  over  the  electrode  area  to  obtain  the  short-circuit  charge,  qsc. 


H.  Integrate  the  deflection,  w,  over  the  face  area;  multiply  by  (0,  and  divide  by  the  face  area 
to  obtain  the  average  face  velocity,  vavg. 

I.  Integrate  the  kinetic  energy  density,  Vi  pv2  =  Vi  p  ofw2,  over  the  ceramic  volume  to 
obtain  the  kinetic  energy,  KE. 

J.  The  blocked  dielectric  energy  density  is  Vi  (£>3  E  3 )s_0.  Use  Eq.  10  to  write  this  energy 
density  and  integrate  over  the  volume  of  the  ceramic  to  obtain  the  blocked  dielectric 
energy,  lfD. 


K.  Determine  the  basic  equivalent-circuit  parameters  - 


m 


d2KE 

dv2 


(18) 
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k 


Cb 

<t> 


d2Usc 

(19) 

d  w2 

avg 

d2UBD 

(20) 

dV2 

dqsc 

3 

(21) 

where  m  and  k  are  the  lumped  mechanical  mass  and  stiffness,  respectively;  Cb  is  the 
blocked  capacitance;  and  0  is  the  electromechanical  conversion  ratio  (current  per  velocity 
or  charge  per  displacement). 


L.  Calculate  the  electromechanical  coupling  factor: 


U°c  -  usc 

K  =  - ~ - 

U°c 

M.  Calculate  the  resonance  frequency: 

2  k 

(0  =  — 
m 


(22) 


(23) 


V.  Solution  Details 

This  section  will  treat  only  the  ceramic  layer.  The  equations  are  complicated  enough  that 
including  the  arbitrary  neutral  plane  location  and  the  substrate  dynamics  will  be  deferred  to 
Appendix  C  and  Appendix  D  as  referenced  at  the  end  of  this  section. 

A.  Find  expression  for  electric  flux  density,  Dj,  in  terms  of  the  electrode  voltage  by 

integrating  Eq.  (10)  over  z.  For  now,  assume  that  the  neutral  plane  is  on  the  inner  surface 
of  the  ceramic. 


t 


V  =  \E3dz  =  -hh3X 


t 

J(5,  +  S2)dz  +  PfctD, 


Solving  for  £>3: 


D,  =  Si 


t 

V  +  hh3l  J(S,  +  S2)dz 


(24) 


(25) 
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B.  Find  the  open-circuit  voltage  by  integrating  Eq.  (24)  over  the  electrode  area  and  setting 
the  surface  charge  equal  to  zero: 

t 

V-A  =  +  S2)dzdA  +  pfctjD^dA  (26) 

AO  A 

The  integral  of  £>3  over  area  (last  term  on  right  side)  is  the  surface  charge.  Set  this  equal  to 
zero  to  obtain  the  open-circuit  voltage: 

t 

voc  =  J('S)  +  S2)dzdA  (27) 

A  0 


C.  Find  Dfc  and  D^c  by  setting  V  =  V°c  and  V=0  respectively  in  Eq.  (25). 


0 


D.  Write  expressions  for  the  strain  energy  density,  U,  for  both  the  open-  and  short-circuit 
electrical  conditions. 


U  =  -5,71  +  -5,7; 

2  2  2 

(30) 

Use  Eqs.  8  and  9  to  replace  T\  and  7V 

uoc 

=  \  {cc°  s!  +  2 ccg  S,  S ,  +  cc"  S2,  -  hhn  D?  (5,  +  S2)} 

(31) 

fjSC 

=  i  {cc,°  S,‘  +  2  cc^  S,  S,  +  cc°  Si  -  hh„  D?  ( S ,  +  .S', ) } 

(32) 

E.  Replace  5|  and  S2  in  the  above  equations  by  their  equivalents  in  terms  of  the  deflection 
function  (Eqs.  16  and  17): 
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Si 

- 

-z 

fa2vvl 

dr1 

= 

-zw' 

(33) 

J 

(idw 

Z  , 

(34) 

^2 

= 

-z 

r  dr 

\  V 

- w 

r 

where  the  primes  indicate  differentiation  with  respect  to  r.  Integration  of  Si  or  S2  with 
respect  to  z  from  0  to  t  is  equivalent  to  replacement  of  z  by  t2/ 2.  Similar  integration  of  either 
quantity  squared  (or  S1S2)  is  equivalent  to  squaring  (or  performing  the  cross  product),  then 
replacing  z2  by  r3/ 3.  (Note:  this  is  where  the  assumption  that  the  neutral  plane  is  at  the 
ceramic-substrate  interface  is  made  explicit;  z  =  0  is  the  location  of  the  neutral  plane  (one 
side  of  the  ceramic)  and  z  =  t  is  the  other  side  of  the  ceramic.)  For  example. 


y°C 

=  hh}'t2  f  (w'  +  w'/r)dA 

2  A  J 

A 

(35) 

Df 

een  hh3l  t  *  ,  ,  v 

= -  (w  +  wlr) 

(36) 

j^OC  _  ^C33  ^31  l 

3  2 

(w#  +  w/r)dA  -  (w'  +  w/r) 

(37) 

A 


F.  Integrate  the  strain  energy  densities  over  the  ceramic  volume  to  obtain  the  total  strain 
energies.  First,  integrate  overz  from  0  to  t: 

t 

J&scdz  = 

0 


\  |jccn  [(w')2  +  (w'lr)1\  +  jcc°2w'w/r  +  hh3]  D*c  (v/  +  w'/r) 
or, 

‘^Usc  dz  =  ^-{/'„[(w')2  +  {w!r)2]  +  /u[2w'w'/r]) 

0 

where  the / parameters  have  been  introduced  to  simplify  the  notation: 


(38) 


(39) 
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3  3 

/'ll  =  CC\\  ~  T®®33  ^31  *  f\\  =  CC\l  ~  "7^33  ^*31 

4  4 
Also,  from  Eqs.  31  and  32, 

t70C  =  Usc  +  ^L(st  +  S2 ) (DjC  -  D°c) 

where,  from  Eqs.  28  and  29, 


D?c  -  D. 


SC 


ee3s3  hhM  t 
2 


+ 


w  / 


(40) 


(41) 


(42) 


so, 


[uocdz  = 

f  Uxiz  -  <K‘‘  («  +  «/r)  f(. 

J 

J  8  A  J  V 

w 


'!  r)d  A 


(43) 


After  integrating  over  the  thickness,  integrate  over  the  cross-sectional  area  to  find  the  total 
energy  density.  Because  only  axisymmetric  deflections  are  being  considered,  the  integration 
over  angle  yields  a  factor  of  2  jt,  consequently, 

Usc  =  J((w')2  +  tylrf)rdr  +  fn  2  jw'w'drj  (44) 

and 

Uoc  =  Usc  +  |  jfr/  +  w'/r)rdr  | *  (45) 


It  is  useful  to  re-write  the  integrals  over  r  in  terms  of  non-dimensional  variables  so  that  the 
integrations  become  functions  only  of  the  shape  of  the  deflection;  actual  dimensions  can  be 
brought  out  into  the  leading  factors.  To  do  this,  define 

£  =  w/w0  ;  T]  s  r/a  (46) 

where  wo  is  the  maximum  deflection  and  a  is  the  disk  radius.  With  these  replacements,  the 
short-circuit  energy  becomes 

V*  =  +  (C'lY'hdn  +  fnljvz'dr,}  (47) 
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with  the  open-circuit  energy  similarly  changed. 

To  further  simplify  the  notation,  it  is  convenient  to  introduce  several  “shape”  integrals  -  that 
is,  integrals  that  depend  only  on  the  shape  of  the  deflection  function: 


(48) 


(49) 


(50) 

(51) 


The  integrals  Ipo  and  Ipi  are  two  terms  in  the  expressions  for  potential  (strain)  energy;  Iq  is 
associated  with  charge;  and  I  a  is  defined  so  that 

A  =  na2IA  (52) 


If  the  electrode  were  continuous  from  center  to  outer  edge,  then  I  a  would  be  one.  For  a 
center-supported  disk  with  a  non-zero  center-post  radius  or  for  an  electrode  that  did  not  cover 
the  entire  disk,  I  a  would  be  less  than  one. 

With  these  shape  integrals, 


rOC 


7Tt3  Wq 
3  a2 


U 


oc 


nt 3  w2 


3  a2 


>1  lQ 

^  hh)\  —  ~ 

(53) 

\f\\  IpO  +  f \ 12  l p\  } 

(54) 

3  / 2  ’ 

0  +  f\2  I pi  2  ^^31  ~  ’ 

(55) 

G.  Integrate  D-iSC  (Eq.  36)  over  the  electrode  area  to  obtain  the  short-circuit  charge,  qsc\ 


qsc  =  ~7iee^hh3{t  j  (w*  +  w'/r)rdr  =  -nw0ee^hh^t  IQ 


(56) 
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Define  a  blocked  capacitance  (that  is,  the  capacitance  with  zero  strain): 


££33  A 


ee^na2  r 

i  A 


t  t 

From  Eqs.  53  and  56, 

=  k cAy0C)1  - 


2  C„ 


and  from  Eqs.  54  and  55, 


Uoc  =  U 


sc 


2  a 2 


+  \Cs(voc)2 


(57) 


(58) 


(59) 


which  shows  explicitly  the  relationship  between  the  pure  electrical  energy  storage  and  the 
two  potential  energy  terms. 

To  complete  the  construction  of  the  equivalent  circuit: 

A.  Integrate  the  deflection,  w,  over  the  face  area  and  divide  by  the  area  to  obtain  the  average 
displacement.  Multiply  this  result  by  co  to  obtain  the  average  face  velocity. 


avg 


! 


2n(o  r  ,  na 

w  rdr  =  — —  CO  wQ  2 


jtridri 


(60) 


=  MWo 


(61) 


where  another  shape  integral  has  been  defined: 

K  =  2^r)dr] 


(62) 


B.  Integrate  the  kinetic  energy  density  over  the  ceramic  volume  to  obtain  the  total  kinetic 
energy,  KE,  in  the  ceramic. 


KE  =  -pv2 
2 


1 


pco 1  w2 


(63) 
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t 

KE  =  —  p<y2  J*J  w2  dz  dA  =  n  pt(02  j  w2  rdr 

AO 

KE  =  pn a2 1  O)2  w]  J  £2  r\dr\  =  pna2 1  (a2  I m 


where 


[  KE 


=  fVvdr, 


(64) 

(65) 


(66) 


C.  Write  the  blocked  dielectric  energy  density  and  integrate  to  find  the  total  blocked 
dielectric  energy.  The  energy  density  is 


U‘°  =  j(A£.) 


5=0 


(67) 


From  Eq.  10  with  Si  and  S2  equal  to  zero: 


A  =  —A 


ee 


33 


(68) 


From  Eq.  25  with  Si  and  S2  equal  to  zero: 


A  =  — V 


(69) 


SO 


£/"  =  \< 


This  expression  is  independent  of  z  and  r  so  the  total  energy  is  simply 

1  eetna 2 1 


UBD  =  -  ee 


1  (V\ 


33 


(J  =  :  ~~33  jA  y  2 


which,  by  Eq.  57,  is 


(70) 


(71) 


U 


BD 


—  CB  V2 
2  B 


(72) 
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D.  Determine  the  equivalent  circuit  parameters.  The  mechanical  mass  is  given  by  Eq.  18: 


m  = 


d2  KE  _  l\  d2  KE 


dvls  !v  d(wvv0)2 


(73) 


With  Eq.  65, 


m  =  2pn a  t 


T2  T 

2 ,  1 A  1  KE 


(74) 


The  mechanical  stiffness  is  given  by  Eq.  19: 


_  d2Usc  _  1 2  d2Usc 

K  —  ^  T 


'w; 


avg 


i:  d 


Wn 


(75) 


With  Eq.  54, 


k  = 


2  nt 


~  2  +  f \ \2  7 P\  )  r 2 

5  a  I 


(76) 


As  an  aside, 


(77) 


The  blocked  capacitance  has  already  been  determined  through  Eq.  57.  Alternatively,  it  is  the 
second  partial  derivative  of  lfD  with  respect  to  voltage.  (Compare  Eq.  71  to  Eq.  57.) 

The  electromechanical  conversion  factor  is  the  ratio  of  electrical  current  under  short-circuit 
conditions  to  the  average  face  velocity.  This  is  identical  to  the  ratio  of  electrical  charge  to 
average  face  displacement: 


V 1  =  [a  Vc 


(78) 


With  Eq.  56, 
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4>  =  /r  ee^  hh3l  t  lnl±  (79) 

1  v 

This  completes  determination  of  the  basic  set  of  equivalent-circuit  parameters.  The  dielectric 
loss  can  be  included  via  the  loss  tangent  and  Cb ■  There  is,  at  present,  no  determination  of 
mechanical  loss.  The  equivalent  circuit9  is  drawn  in  Fig.  5  using  an  explicit  electrical-to- 
mechanical  transformer  and  in  Fig.  6  by  transforming  all  elements  to  the  electrical  side. 


Fig.  5.  Electromechanical  equivalent  circuit  for  flexural  disk  transducer  with  mechanical 
quantities  to  the  right  of  the  transformer  and  electrical  quantities  to  the  left  of  the 
transformer. 


Fig.  6.  Equivalent  circuit  with  all  quantities  expressed  in  equivalent  electrical  terms. 


(See  Eqs.  59  and  77.)  With  Eqs.  54  and  55, 
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K  = 


-  ee3?3  hh]x  — 

2  IA 


f\\  I pO  +  f 12  I, 


/»■ 


+  -4K,  — 

2  /, 


(81) 


F.  Calculate  the  resonance  frequency: 

0)]ir  =  A:  /  m 


(82) 


In  this  version,  radiation  loading  is  considered  by  including  the  approximate  induced  water 
mass.  The  water  mass  is  calculated  from  the  mass  appropriate  to  a  uniform,  baffled  piston  of 
the  same  effective  area  as  the  actual  piston.  This  radiation  mass  is 


m 


rad 


=  P* 


na 


f%a} 
3  n 


(83) 


The  water-loaded  resonance  frequency  is  then 


(0 


2 

water 


k 

m  +  mrad 


(84) 


The  infinite  baffle  condition  was  chosen  since  it  is  a  reasonable  approximation  in  two 
common  situations.  The  obvious  case  is  if  the  transducer  is  mounted  in  a  larger,  nearly  flat 
structure.  Less  obvious  is  the  case  in  which  the  transducer  comprises  two  flexural  disk 
elements  mounted  back-to-back  and  driven  to  expand  and  contract  in  phase.  In  this  case,  the 
pair  of  sources  can  be  replaced  by  a  single  source  in  an  infinite  baffle  since  the  pair  of 
sources  represents  the  image-source  solution  to  the  baffled-source  problem.  Even  if  the 
element  is  not  baffled,  the  error  introduced  is  relatively  small. 

Treatment  of  an  arbitrary  neutral  plane  and  incorporation  of  the  dynamics  associated  with  the 
substrate  are  straightforward  modifications  of  the  above  analysis.  The  details  are  included  in 
Appendix  C,  Appendix  D,  and  Appendix  E. 

A  preliminary  treatment  of  stress  in  the  transducer  element  is  outlined  in  Appendix  F  and  the 
results  of  the  analysis  in  this  report  are  compared  with  the  results  given  in  Woollett’s  report1 
in  Appendix  G.  Appendix  H  contains  the  MathCad  worksheet  listings  and  Appendix  I 
contains  the  materials  properties  files. 
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Apppendix  A:  Exact  Deflection  Functions 

For  thin-plate  theory,  the  general  solution  for  the  transverse  displacement,  w,  as  a  function  of 
radius,  r,  and  angle,  6,  is 

w(r,6)  =  cos  (nd)[AniJn{kr)  +  AnlYn{kr )  +  AniIn(kr)  +  An4Kn(kr )]  (Al) 

where 

*<  =  12(‘-CTV;  (A2) 

Et2  v  ' 

and  E  is  the  material’s  Young’s  modulus  and  a  is  the  Poisson’s  ratio.  J„  and  Y„  are  the  ordinary 
Bessel  functions  of  the  first  and  second  kind,  while  I„  and  K„  are  the  modified  Bessel  functions. 
Here,  only  the  axially  symmetric  solutions  will  be  considered  ( n  =  0): 

w(r)  =  AxJ0{kr)  +  A2  Y0  (kr)  +  A3I0(kr)  +  A4K0(kr)  (A3) 

To  organize  the  boundary  conditions,  a  number  of  auxiliary  definitions  will  be  introduced. 

These  include  displacement  functions, 

D\  =  Jq  >  D2  —  Y0  ;  D3  =  I0  ;  D4  =  K0  (A4) 


slope  functions, 

5,  =  -J,  ;  S2  =  ~  Y{  ; 

S,  =  /,  ;  S,  =  -K, 

(A5) 

radial  moment  functions', 

js: 

n 

> 

M *  '  \rY' 

(A6) 

v  * 

1 

II 

cn 

M •  '  \rK' 

(A7) 

and  radial  (Kelvin-Kirchhoff)  shear  functions, 

v  =  -  /  •  V  =  -  Y 

r\  J\  »  vi  L\  ’ 

V3  =  -I,  ;  V4  =  Kt 

(A8) 

A-  1 


Flexural-Disk  Transducer  Model 


T.  B.  Gabrielson 


These  functions  are  combined  to  form  the  appropriate  boundary  conditions.  For  example,  the 
condition  for  zero  displacement  at  r  =  a  would  be  written  as  follows: 


J^A,D,(ka)  =  0  (A9) 


The  frequency  constant  is  found  by  searching  for  the  zeroes  of  the  determinant  of  the  matrix 
formed  by  the  boundary  condition  functional  forms.  For  edge-supported  plates,  there  will  be  two 
conditions  on  the  edge  and  the  matrix  will  be  2x2;  for  plates  supported  by  an  inner  support  post, 
there  will  be  two  conditions  at  the  support  and  two  conditions  at  the  outer  edge  so  the  matrix  will 
be  4x4.  If  a  is  the  outer  radius,  then  define  a  frequency  constant.  A2,  equal  to  ( ka )2.  For  a 
clamped  edge,  the  deflection  ( D )  and  the  slope  (S)  will  be  zero;  for  a  simply-supported  edge,  the 
deflection  ( D )  and  the  radial  moment  (M)  will  be  zero;  for  a  free  edge,  the  radial  moment  (M) 
and  the  shear  (V)  will  be  zero. 


For  example,  the  boundary-condition  matrix  for  a  circular  plate,  simply-supported  at  the  edge  is 


£>,  (A)  £>3  (A)' 
A/,  (A)  M3(A) 


(A10) 


The  functions  with  subscripts  2  and  4  are  associated  with  the  Bessel  functions,  Y,  and  K,  which 
are  infinite  for  zero  argument;  consequently,  these  functions  are  discarded  in  the  solution  for 
plates  with  continuity  at  r  =  0. 

As  another  example,  the  boundary-condition  matrix  for  a  circular  plate,  centrally  clamped  at  r  = 
b  and  free  at  the  outer  edge  is 


D ,  (aA) 

D2  (a A) 

Z)3  (aA) 

D,  (aA)' 

5,  (aA) 

S2  (aA) 

S2  (aA) 

St  (aA) 

M,(A) 

mM 

M,(A) 

M.(A) 

.  r,(A) 

nw 

V .  W  . 

(All) 


where  a  =  b/a.  In  matrix  form,  the  actual  boundary  conditions  would  be  written 

F-A  =  0 


where  A  would  either  be  2x1  or  4x1 : 


A 


edge 


A 


A 


A 


center 


A, 

A 

A 


A 


(A  12) 


(A13) 
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In  any  case,  the  frequency  constants  are  the  values  of  A  for  which 

1^1  =  0  (A  14) 


The  frequency  constants,  Am,  are  independent  of  the  absolute  disk  dimensions  and  the  disk 
material  except  for  a  modest  dependence  on  Poisson’s  ratio  (through  the  M  functions).  There  are 
an  infinite  number  of  axially  symmetric  resonances  but,  here,  only  the  first  is  important  so  the 
search  for  the  roots  of  Eq.  A 14  would  be  restricted  to  the  smallest  value  of  A. 

Once  the  frequency  constant  for  the  lowest  mode  is  found,  the  deflection  function  can  be 
determined.  For  the  edge-supported  cases,  there  are  two  arbitrary  constants  (A  \  and  A3  from  Eq. 
A 13).  At  resonance,  the  two  equations  resulting  from  the  matrix  multiplication  in  Eq.  A12  are 
not  independent  -  either  represents  the  solution.  One  of  those  equations  determines  the  ratio  of 
A 1  and  A3.  Normalizing  the  equation  so  that  the  maximum  deflection  is  one  completely 
determines  the  constants. 

For  the  centrally  supported  case,  there  are  four  arbitrary  constants.  Determination  of  the 
constants  proceeds  as  follows:  (1)  Set  ^4 1  to  one.  (2)  Add  the  first  two  equations  that  result  from 
the  matrix  multiplication  in  Eq.  A12  together.  (3)  Form  the  3x3  system  of  equations  from  step  2 
and  the  remaining  two  equations  with  the  constants  A2-A4  on  one  side.  (4)  Solve  this 
inhomogeneous  equation  system  for  A2-A4.  (5)  Normalize  the  result  so  that  the  maximum 
deflection  is  one. 


A  -  3 
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Appendix  B.  Rayleigh-Ritz  Approximation 

Exact  solution  for  the  deflection  function  in  the  thin-plate  approximation  is  straightforward  as 
long  as  the  disk  is  uniform  in  composition  and  thickness.  In  practice,  however,  the  ceramic  in 
most  flexural  disk  transducers  does  not  completely  cover  the  substrate  and  this  can  have 
important  consequences  especially  with  regard  to  inner-edge  stresses  and  electric-field 
breakdown.  To  anticipate  the  need  to  treat  flexural-disk  structures  in  which  the  ceramic  only 
partially  covers  the  substrate,  an  approximate  energy-based  technique  for  finding  the  deflection 
function  is  introduced  in  this  second  release.  This  solution  is  one  of  the  Rayleigh-Ritz 
techniques.  In  the  polynomial  form  of  the  Rayleigh-Ritz  technique,  the  deflection  function  is 
assumed  to  have  the  form  of  a  polynomial  (normally  with  degree  between  two  and  ten).  The 
resonance  frequency  is  computed  from  the  strain  energy  and  the  kinetic  energy  and  then  the 
resonance  is  minimized  with  respect  to  the  coefficients  of  the  polynomial.  In  this  manner, 
several  solutions  are  produced.  Each  solution  is  an  approximation  for  a  mode;  the  lower  the 
mode,  the  better  the  approximation.  Generally,  the  fundamental  mode  is  approximated  quite 
closely  while  the  highest  modes  are  approximated  poorly. 

The  starting  point  for  the  Rayleigh-Ritz  technique  is  construction  of  a  function  that  satisfies  the 
boundary  conditions  and  that  has  a  number  of  adjustable  parameters.  Then  both  the  potential 
energy  and  the  kinetic  energy  are  expressed  in  terms  of  this  function.  The  resonance  frequency 
is  obtained  by  setting  the  kinetic  energy  equal  to  the  potential  energy  and  solving  for  frequency. 
Minimizing  the  frequency  with  respect  to  all  of  the  adjustable  parameters  produces  an 
approximate  solution  for  both  the  resonance  frequencies  and  the  mode  shapes.  One  of  the 
simplest  implementations  for  the  Rayleigh-Ritz  function  is  as  a  polynomial. 


Application  to  transverse  vibrations  of  a  circular  membrane 

To  illustrate  the  polynomial  Rayleigh-Ritz  solution,  we  will  start  with  a  simpler  problem:  axially 
symmetric  transverse  vibrations  of  a  circular  membrane.  After  this  introductory  illustration,  the 
technique  will  be  applied  to  flexural  vibrations  of  a  thin  circular  plate. 

First,  write  the  displacement  of  the  membrane  as  a  polynomial  in  the  radial  coordinate.  Here,  we 
will  use  normalized  forms.  The  variable,  £  is  the  displacement  divided  by  the  peak 
displacement  (at  the  center  of  the  membrane).  The  variable,  77,  is  the  radial  coordinate  divided 
by  the  radius  of  the  membrane. 


£(*?)  =  a0  +  a,  7]  +  a2r) 2  +  a3? f  +... 

The  boundary  conditions  for  the  membrane  are:  (1)  displacement  equal  to  zero  at  the 
circumference,  and  (2)  slope  equal  to  zero  at  the  center: 
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(Actually,  the  condition  on  the  slope  at  the  center  is  unnecessary.  The  accuracy  of  higher-order 
modes  is  slightly  better  using  both  conditions.) 

The  slope  is 

— —  =  <2|  +  2  a2  r)  +  3  a3T)2  +... 
dri 


so  we  set 


and 


dr) 


=  a,  =  0 


n=  o 


£(l)  =  <20  +  +  <32  +  Oj  +  ...  =  0 

•  •  Gq  —  —  @2  —  ~~  @4 

Therefore,  the  normalized  deflection  can  be  written  as  follows: 

£(»?)  =  -  a2  +  a2r\ i2  -  a3  +  a3rf 


or 


-  aA  +  a4r)  -  ... 

=  (V2  ~  l)  +  a,  (t?3  -  0 

+  °4  (r?4  ~  l)  +  "• 

Notice  that,  in  the  second  form,  each  term  satisfies  the  boundary  conditions. 
For  convenience,  we  change  the  indexing  to  1  ...  N  and  replace  am+ 1  by  bm: 


f(n)  =  £  “  ') 


(-♦dm* 


The  kinetic  energy  of  the  membrane  is 
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or,  in  terms  of  the  normalized  displacement, 


KE  =  o2  w2^  pna2 1  \^2  r]dr\ 

o 

The  potential  energy  is  equal  to  the  work  against  the  membrane  tensile  stress,  <7.  For  an  element, 
dr,  of  the  membrane,  the  unstretched  length  is  dr  while  the  stretched  length  is  the  square  root  of 
dr2  plus  dw2.  The  elongation  is  the  difference  between  these  two  lengths  which,  for  small 
deflection,  is  approximately 


^  dw^ 
\dr  ) 


dr 


The  force  against  which  this  elongation  takes  place  is  the  product  of  the  tensile  stress  and  the 
area  of  the  edge  of  the  annular  element  of  the  membrane  so 


d  PE  =  Inrta  — 
2 


^  dw^ 

ydr, 


dr 


The  total  potential  energy  (at  maximum  deflection)  is  then 


PE  =  h \Kknat\ 
0 

If  we  define  U  and  T as  follows: 

u  ,  j 


Kdn, 


r]  dr] 


1 

f 

(d£\ 

J 

0 

A 

r]dr)  ;  T  =  \^r]dr] 

0 

we  can  equate  the  kinetic  and  potential  energies  and  solve  for  the  frequency: 


2  a 
or  =  — ■ 


pa 


U_ 

T 


We  can  also  define  a  normalized  resonance  frequency  for  membrane: 

U 


a  = 


B  -  3 
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To  find  the  minima  with  respect  to  bm: 

d^_  =  T(duldb.)  ~  u(dTldb.)  =  » 
db.  T1 


Set  the  numerator  (divided  by  7)  equal  to  zero: 


which  is  the  same  as 


dU  U  dT 

dbm  T  dbm 


dU  2  dT 

-  -  a  - 

db  db_ 

m  m 


0 


This  has  the  form  of  a  matrix  eigenvalue  problem.  In  matrix  terms,  if  B  is  a  column  vector  of  the 
coefficients,  bm,  then 


T  =  B'  T  B 


and 


—  =  [0  0  ...  1  0]T  B  +  B'r  T  [0  0  ...  1  Of 

dK 

where  the  row  vectors  are  zero  except  for  a  one  in  the  wth  position.  The  column  vector  of  the 
derivatives  is  then 

=  TB  +  (B'rT)  =  (T  +  T,r)B 

Therefore, 


dT 


T  =  T  +  fr 


and  similarly  for  U. 
Since, 


dU  dT 
dbm  '  db 

m  m 


^  bn  x  functions  of  n,m 
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there  are  N equations  each  involving  the  N  values  of  bm\ 

(NxN)  ( jVxl ) 

U  B  -  a2  T  •  B 


0 


where 

and 


U ,  T  =  NxN  square  matrices 


Nxl 


The  equations  are  homogeneous  so  solutions  only  exist  for  specific  values  (eigenvalues)  of  a2: 

(U  -  a2  T)  •  B  =  0 

Each  a2  has  an  associated  B  (the  eigenvector  or  column  vector  of  bm’s). 

So  a2  could  be  written  as  a  diagonal  matrix: 


[«•]  - 


a; 


a 


22 


a 


NN 


Since, 


U  -  a2  T  =  0 
a2T  =  U 


[a2  ]  =  diagonalization  of  T~‘  •  U 


If  we  calculate  T'1  •  U  ,  then  a  standard  matrix  eigenvalue  solver  can  be  used  to  find  both  the 
eigenvalues  and  the  eigenvectors. 

Continuing  the  example  of  the  membrane  in  tension,  U  and  T  are  found  as  follows.  The 
integrand  of  the  kinetic  energy  is 

t‘n  =  JX  b,b,(n-'  -  i)  (, r'  -1)1 

m  n 

=  D nm+n+ 3  -  r)m+1  -  T)n+2  -  7]] 
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Integrate  to  find  T: 


T 


E  E  KK 


1 

m  +  n  +  4 


1 

m  +  3 


then  differentiate  with  respect  to  bm : 


§ -  22>- 

dbm  t) 


1 


1 


/m  +  m  +  4  /m  +  3 


1 

m+3 


1 

m  +  3 


2 


1 

2 


The  factors  in  square  brackets  are  the  matrix  elements  of  T  .  Notice  that  T  is  symmetric 
hence,  when  added  to  its  transpose  merely  produces  the  leading  factor  of  two;  therefore,  the 
elements  of  T  are  two  times  the  quantities  in  square  brackets. 

The  integrand  of  the  potential  energy  is 


(g}1 


V  =  E  X  bmb„(m  +  1)  {n  +  1  )r] 


m+n+l 


so  integrate  to  find  U: 


U 


1 

f 

(d£ 

J 

0  1 

r]dr ? 


E  E 


(m  +  1)  (m  +  1) 
mj  +  m  +  2 


and  differentiate, 


dU 

dbm 


(m  +  1)  (m  +  1) 
m  +  m  +  2 


where  the  factors  in  square  brackets  are  the  matrix  elements  of  U  . 


B  -  6 
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Application  to  transverse  vibrations  of  a  center-supported  circular  disk 

For  the  centrally  supported  annular  disk,  the  boundary  conditions  to  be  satisfied  are  for  the 
displacement  and  slope  at  the  inner  edge  to  be  zero.  In  energy-minimization  methods,  it  is  not 
necessary  to  satisfy  explicitly  conditions  on  moment  or  shear  (at  hinged  or  free  boundaries,  for 
example)81.  The  process  of  minimization  automatically  satisfies  these  conditions  on  derivatives 
higher  than  the  first.  Note  that  in  this  section,  b  refers  to  the  inner  radius  (the  radius  of  the 
clamp);  it  is  not  a  polynomial  coefficient. 

The  polynomial  solution  that  satisfies  the  clamped  inner-edge  boundary  conditions  is 

£  =  5X[y  +  1  -  {m  +  l)bmv  +  rnbm+ '] 

The  first  and  second  derivatives  are 

=  5>.  (»+ 0[>i‘  -*■] 

We  have  already  treated  the  normalization  constants  in  the  first  release  report  so  we  need  only 
construct  the  mode  shape  using  Rayleigh-Ritz.  Once  we  have  the  normalized  shape,  we  can  use 
the  shape  integrals  previously  derived  to  obtain  the  equivalent-circuit  parameters  and  resonance 
frequency.  Consequently,  to  find  the  shape,  we  can  ignore  all  leading  factors  and  assume  an 
effective  modulus  and  Poisson’s  ratio  for  the  composite  disk. 

Since 

Tx 

t2 

the  strain  energy  density  is 

Usr  =£*,»;  +  jS.  T, 

=  5T^[S'  +  +  2v^J 

81  R.  Weinstock,  Calculus  of  Variations,  Dover,  NY,  1974,  §7.7. 


-  +  «-%] 

=  +  ^.1 
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Although  this  form  is  useable,  it  is  commonly  rearranged  as  follows: 

1  E 


U, 


ST 


=  2T7\-{s'  +s')2  "  2(1~v)s'sk 


Substituting  for  the  strains, 


Usr  — 


2  1  —  v2 
z2  E 


f  .  w'Y 

-  +  7 


2(1-,) 


0  t 

w  w 


2  1  —  v2 


\d_ 

r  dr 


dw 

dr 


\2 


)) 


V  ’  r  dr1  dr 


The  integration  over  the  volume  of  the  composite  disk  produces  several  leading  factors  that  can 
be  ignored.  The  integration  over  the  radial  coordinate  is  necessary  to  retain  but  all  the  shape 
information  is  contained  in  the  normalized  forms.  Therefore,  let 


u  =  Lm  -  /„ 


where 


I  sum 


fl 

r  d 

(dt\  1 

in 

dn)\ 

dr) 


and 


=  2J 


L 


dr)2  dr) 


dr) 


(The  notation  of  this  last  integral  is  consistent  with  that  in  the  first  report.  ISUm  is  equal  to  the  sum 
of/po  and/pi.) 


Each  integral  produces  a  matrix.  The  algebra  is  straightforward  although,  in  some  cases,  tedious. 
To  illustrate,  consider  the  integrand  of  Isum. 


dr) 


d£ 

V~r 


\  “|2 


=  X5X  a„  (m  +  l)(/i  +  l)> 


(m  +  l)(n  +  \)t)m+n-'  -  bm  (n  +  \)i)n-'  -  bn  (m  +  l)r)m-'  + 


V 


B  -  8 


Flexural-Disk  Transducer  Model 


T.  B.  Gabrielson 


The  integration  produces  the  double  sum  of  the  product  of  am,  a„,  and  the  matrix  elements  of  the 
desired  matrix: 


(m  + 1)  (n  + 1) 
m  +  n 


(m  +  l)(/j  +  l)» 

»  +  l  bm 
n 


(l -6n+") 


(*-»') 


— *■ 
m 


(>-»■) 


bm+n  In  b 


In  general,  to  produce  the  required  matrices,  we  can  take  the  complete  factor  that  multiplies  the 
product  of  am  and  a„  in  the  double  summation,  divide  it  by  two  and  add  it  to  its  transpose.  If  the 
factor  is  symmetric  in  m  and  n,  though,  this  step  is  redundant. 

The  result  for  Ipi, 


=  X5X  an  2n(m  +  l)(n  +  \) 

is  not  symmetric  in  m  and  n  so  we  would  generate  one  matrix  from  half  of  the  factor  shown,  then 
add  it  to  its  transpose  to  find  Up|.  The  complete  matrix  for  potential  energy  for  the  disk  is  then 

u  =  U_  -  (l-v)U„ 

The  matrix  for  kinetic  energy  is  found  in  similar  fashion: 


1  -bm+n  bM(\-bn) 

m  +  n  n 


T(m,n) 

+  nb 


1  -  b" 


m  +  n  +  4 
■  1  -  bm+i 


-  (n  +  1)6' 


1  -  b‘ 


i+4 


m  +  4 
1  -  6"+4 


-  (m  +  1)6' 

m+3  7  n+ 4 


+  (m  +  1)(«  +  l)6m+n^ — -  ( 2mn  +  m  +  n)bm+n+'  - — 


+  mb 


n+  3 


4 

+  mn  b" 


-+.  l-^n+3  .  -  b1 


Once  U  and  T  have  been  found,  the  matrix  eigenvalue/eigenfunction  solver  is  used  to  find  the 
fundamental  mode  frequency  and  mode  shape.  Matrix  eigenvalue  solvers  do  not  usually  return 
the  results  in  order  from  lowest  to  highest  mode;  the  smallest  eigenvalue  corresponds  to  the 
fundamental  mode.  The  MathCad  rsort  function  is  used  so  that  both  the  eigenvalue  matrix  and 
the  eigenvector  matrix  can  be  sorted  at  the  same  time. 

Once  the  mode  shape  has  been  found  (that  is,  the  polynomial  coefficients  have  been  determined 
from  the  eigenvector  values),  the  deflection  can  be  normalized  so  that  the  maximum  deflection  is 
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one  for  consistency  with  the  first  release.  Then  the  other  shape  integrals  can  be  calculated.  If  the 
normalized  polynomial  coefficient  matrix  is  C,  then, 


K  =  =  C'TC 


I,,.,  -  C"U„,C 

=  C"  u„  c 

Ip  0  I psum  I p\ 


As  before,  the  shape  integral  associated  with  the  effective  area  is 

/„  =  2jt]dn  =  i  -  b2 


The  shape  integral  associated  with  surface  charge  is 


1 

/ 


v  + 


r\dr\  =  C’Q 


where 


Q(m)  =  (w+l)(l-6m) 
and  the  shape  integral  associated  with  velocity  is 

Iv  =  2  j^dn  =  c,rv 


where 


V(m)  =  2 


1  —  Am  +  3  1  —  h* 

— — -  -  (rn  +  l)bm  — -  +  mb " 

m+ 3  V  ;  3 


+ll~b2 


For  calculation  of  the  deflection  corresponding  to  maximum  strain  (see  Section  III),  the 
following  definition  is  useful: 


SM  (m)  =  m(m  +  1  )bm-' 


so  that 


£'(b)  =  C'SM 
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The  matrix  forms  permit  a  cleaner  program  since  U,  T,  Q,  and  V  can  be  precomputed.  Then  the 
shape  integrals  are  simple  matrix  multiplications  with  the  coefficient  matrix,  C. 


Application  to  transverse  vibrations  of  a  edge-supported  circular  disk 

For  the  edge-supported  disk  (with  simply-supported  outer  edge),  the  necessary  boundary 
conditions  are  that  the  deflection  be  zero  at  the  outer  edge  and  that  the  slope  be  zero  at  the 
center: 


4(1)  =  o  ;  «'(o)  =  o 

The  normalized  deflection  function  and  its  derivatives  are  then 

4  -  X  k  (r"  -  i) 

4'  =  +  i)i" 

4'  =  +  !)!,-' 


The  required  matrices  are 


and 


u  (m  n)  _  (m  +  O2  (n  + 


(m , «) 


n  (m  +  l)  (n  +  1) 
(m  +  n) 


T  (m,n)  = 


u,.  =  +  u-;, 


u  =  u_  -  (l-v)upl 

1  1  1 


m  +  n  +  4  m  +  3  n  +  3 


2 


Q(m)  =  m  +  1 


v(m)  = 


m  +  3 


SM(w)  =  m  [m  +  1)  & 


m  - ! 
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The  formulation  in  this  release  is  sufficiently  general  that  the  edge-supported  structure  can  be 
modeled  either  with  or  without  a  substrate.  (Set  the  thickness  of  the  substrate  to  zero  in  the  latter 
case.)  The  edge-supported  structure  is  assumed  to  be  two  identical  ceramic  disks  laminated  to 
each  other  or  on  opposite  sides  of  a  substrate  disk. 
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Appendix  C.  Arbitrary  Location  of  Neutral  Plane 

While  there  is  some  logic  in  locating  the  neutral  plane  at  the  interface  between  the  ceramic  and 
the  substrate,  this  specific  choice  limits  the  design  excessively.  A  module  was  added  to  compute 
the  location  of  the  neutral  plane  for  arbitrary  thickness  of  ceramic  and  substrate  (within  the  limits 
of  the  thin-plate  approximation). 

The  neutral  plane  can  be  found  by  minimizing  the  strain  energy  in  the  two-layer  system.  If  the  z- 
coordinate  is  defined  so  that  z  =  0  is  the  location  of  the  neutral  plane,  then  the  strain  in  the 
system  has  the  following  relationship  to  the  deflection,  w: 

d2w 

z — 7-  =  -zw 
dr2 

z  dw  _  zw 
r  dr  r 


^ - 

a 

b*a 

r  1 

t, 

SUBSTRATE  1 

Pit 
'  if 

■  - 

— U- 

CENTER 

CERAMIC  A 

tc 

SUPPORT  POST 

| 

Fig.  Cl.  Center-supported  flexural  disk  structure.  In  this  view,  the  water  would  be  on  the  upper 
face  and  the  lower  face  would  contact  a  gas  back-volume.  For  many  practical  transducers,  the 
structure  would  be  doubled  so  as  to  have  two  structures  back-to-back.  The  edges  are  nominally 
free  but  would,  in  practice,  be  sealed  with  boots  or  bellows. 

Let  the  ceramic  be  below  the  substrate  and  the  ceramic  thickness  and  substrate  thickness  be  tc 
and  ts  respectively  (see  Fig.  Cl).  If  the  neutral  plane  is  at  z  =  0  and  we  assume  that  the  neutral 
plane  is  in  the  substrate  a  distance  zq  from  the  interface  (see  Fig.  C2),  then  the  ceramic  extends 
from  z  =  -  tc-zo  to  z  =  -  zo,  and  the  substrate  extends  from  z  =  -  zo  to  z  =  ts  -  zq. 
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z  -  ts  -  z  o 

z  =  0 
z  =  -z0 

z  =  -2o  -  *c 

Fig.  C2.  Relationship  of  neutral  plane  (dash-dot  line)  to  substrate/ceramic  interface.  The  origin  of 
the  z-coordinate  system  is  taken  to  be  at  the  neutral  plane. 

The  strain  energy  density  is 

0„  =  \s"t  =  is,7;  +  is,r2 

since  all  stresses  other  than  T\  and  Ti  are  zero. 

The  strain  energy  is  the  integral  of  the  strain  energy  density  over  the  entire  structure.  Integration 
over  r  and  9  produce  factors  that  do  not  depend  on  the  location  of  the  neutral  plane,  so  these 
integrations  do  not  need  to  be  carried  out  explicitly.  Since  the  faces  of  the  ceramic  are 
electroded,  the  electric  flux  density,  D,  does  not  depend  on  z  so  we  will  use  the  basic  equation: 


which  expands  to 

Tx 

•  t2 
0 


T 

= 

cDS 

cRSx 

+ 

c°S2 

+ 

C13 

+ 

cn  S2 

+ 

c°Sx 

+ 

<&s2 

+ 

C3°3  S3 

The  last  equation  can  be  used  to  eliminate  S3  from  the  first  two  to  produce: 


f 

2  > 

f 

2  X 

zz 

^11 

_  £n 

st 

+ 

C\2 

_  £ii 

V 

C33  ) 

V 

C33  j 

For  convenience,  redefine  the  quantities  in  parentheses  as  the  effective  constants,  ccD: 


Tx  =  ccxx  S,  +  ccx2  S2 

Also, 

T2  =  ccX2  S,  +  cc°  S2 
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(Note:  The  cc  coefficients  can  also  be  obtained  by  the  matrix  relations  discussed  previously. 
While  the  forms  are  different  - 


cc,. 


cc,- 


the  values  are  identical. 


For  an  isotropic  material,  the  Young’s  modulus,  E,  is  the  reciprocal  of  s  1 1  and  the  Poisson’s  ratio, 
v,  is  negative  the  ratio  of  S12  to  sn.  This  produces  the  usual  forms  for  thin  plates: 


Tx 

t2 


+ 


+ 


It  is  useful  to  remember  that,  while  T3  is  assumed  to  be  zero,  S3  is,  in  general,  not  zero.) 

The  strain  energy  density  is  then 

UST  =  ^(s,2  +  S]  +  2vS,S2) 

(We’ll  drop  the  superscript,  D,  on  cc  so  that  we  can  use  the  superscript  to  differentiate  between 
the  ceramic  and  the  substrate.) 

Actually,  the  strain  energy  also  contains  terms  with  the  product  of  S  and  D.  We  will  ignore  these 
terms  with  respect  to  locating  the  neutral  plane.  Because  the  stiffness  of  a  piezoelectric  material 
depends  on  the  electrical  boundary  conditions,  the  location  of  the  neutral  plane  also  must  depend 
on  these  conditions.  For  the  present,  this  is  presumed  to  be  a  small  effect;  however,  if  the 
coupling  factor  is  large,  it  may  be  necessary  to  modify  this  assumption.  We  can  test  the 
assumption  by  calculating  the  neutral-plane  location  using  both  cc°  and  ccE  coefficients. 

Each  of  the  terms  in  the  strain  energy  density  equation  is  the  product  of  two  strains. 
Consequently,  each  term  has  as  a  factor  z2.  This  z  factor,  cci  1,  and  v  are  the  only  quantities  that 
change  with  z  (cc  and  v  because  the  material  changes  in  the  thickness  direction).  We  will 
further  assume  that  the  Poisson’s  ratios  of  the  two  materials  are  not  radically  different.  We  will 
retain  the  dependence  of  cc  on  v  but  ignore  the  change  in  v  as  a  factor  of  S1S2. 
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Integration  over  r  and  6  produce  factors  that  have  no  dependence  on  the  location  of  the  neutral 
plane;  these  factors  can  be  ignored  in  the  energy  minimization.  Integration  over  z  produces  a 
factor  of  z3  times  cci  1.  Integration  over  the  ceramic  extends  from  z  =  -  tc  -  zq  to  z  =  -  zo  to 
produce 


«*r [(-*.)’  -  k  -hf 


while  integration  over  the  substrate  extends  from  z  =  -  zo  to  z  =  ts  -  zo,  producing 

ccu  [(*,  -  zo)3  -  (-z0)3] 

If  we  add  these  two  expressions  together  (for  the  total  strain  energy),  differentiate  with  respect  to 
zo,  set  equal  to  zero,  and  solve  for  zo,  we  obtain 


+  <<»,] 


This  gives  the  location  of  the  neutral  plane  in  terms  of  the  material  properties  and  thicknesses  for 
the  ceramic  and  substrate  layers.  If  zo  is  negative,  then  the  neutral  plane  is  in  the  ceramic  and  a 
hydrostatic  load  would  put  part  of  the  ceramic  in  tension  unless  prestressed. 
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Appendix  D.  Accounting  for  the  Substrate 

In  the  first  release,  the  neutral  plane  was  assumed  to  be  at  the  interface  between  the  ceramic  and 
the  substrate  and  the  substrate  properties  were  taken  to  be  identical  to  the  ceramic  properties. 
This  meant  that  only  the  ceramic  needed  to  be  analyzed.  In  subsequent  releases,  arbitrary 
properties  and  thickness  were  permitted  for  the  substrate.  Consequently,  the  calculations  in  the 
ceramic  were  modified  (since  the  neutral  plane  is  no  longer,  in  general,  at  the  inside  face  of  the 
ceramic  and  calculations  in  the  substrate  were  added. 

The  primary  difference  in  the  mechanics  of  calculation  is  that  the  limits  of  integration  in  the  z- 
direction  are  no  longer  zero  and  tc.  Integrations  over  thickness  in  the  ceramic  go  from  -zo~tc  to 
-  zo  and  integrations  over  thickness  in  the  substrate  go  from  -  zo  to  ts  -  zo.  Before,  integrations  of 
S i  or  S2  over  the  ceramic  produced  a  factor  of  tc  /  2;  now  the  same  integrations  produce  a  factor 
of 


1  + 


2  zn 


where  we  have  defined  a  constant,  ac\,  equal  to  the  factor  in  square  brackets.  Similar  integration 
over  the  substrate  produces  a  factor  of 


1  +  ^ 


=  —a. 


If  the  integrations  are  of  the  squares  of  strains,  then  a  factor  of 


1  + 

L 


is  produced  over  the  ceramic  and  a  factor  of 

3  zn 


+ 


9  zn 


*e  2 


1  - 


9z„ 


_  ^ 


s2 


is  produced  over  the  substrate.  (Notice  the  minus  sign  in  this  factor.) 

Several  quantities  are  unaffected  by  these  changes.  The  blocked  electrical  capacitance  is  still 


C,  - 
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the  effective  area  is 


and  the  average  face  velocity  is 


A  =  na 2 1. 


V  =  Wn  CO  — 

avg  0  j 

1 A 

The  open-circuit  voltage  changes  to 

V™  =  -^MhxLLaA 

The  potential  energy  in  the  ceramic  for  short-circuit  electrical  conditions  becomes 


Usc  = 

cer  2 

a  3 


'c2  [/ll  IP0  + 


fn  Ip\  ] 


where 


f\\  —  CCU  A  ^33  ^^31 

4 


/  V 


y°c2  j 


,  fl2  ccx 2  ££33  hh} 


(  „  V 


v  °c2  J 


and  the  potential  energy  in  the  ceramic  for  open-circuit  electrical  conditions  is 


Uoc  =  Usc  + 

cer  cer 


The  potential  energy  in  the  substrate  is 

U  -  2^. a  f  csub  j  +  csub  /  1 

u sub  Q2  2  Usl  l  11  1p°  T  0,2  lP'\ 

where 


„sub 


“'sub 


1  - 


>  C\2 


v  csub 

vsub  cll 
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The  total  charge  under  short-circuit  conditions  is 

q  =  w0  ££33  hhil  n  tc  <2fl  I A 


and  the  kinetic  energy  is 


KE  =  wl  no1  o2 [pc  tc  +  psts ] 

The  elements  of  the  equivalent  circuit  are  defined  in  the  same  way  but  have  somewhat  different 
values.  The  equivalent  mass  is 


m 


d2KE 


/  /, 


=  ^  =  2  na2^M-[pctc  +  Plt,] 

y  *ayg  *v 


The  equivalent  stiffness  is 


d2Usc 


k  = 


3a2  I2  ^  *c  a<:2  Ip0  +  f'2  ^  +  **  °s2 f-C”*  Ip0  +  ^  I p'  ]} 


3  w‘ 


avg 


where  the  potential  energy  is  the  sum  of  the  energy  in  the  ceramic  and  the  energy  in  the 
substrate.  The  transduction  factor  is 


<t> 


dq 


sc 


3  w. 


avg 


7t  ££33  /t/iji  tc  acl 


IqIa 

K 


The  coupling  factor  is  most  easily  calculated  using  the  following  form: 


K 


2 


U°C  _  ySC 


To  estimate  an  upper  limit  to  the  transmitting  voltage  response,  we  assume  operation  at 
resonance  and  assume  that  the  mechanical  loss  in  the  transducer  is  much  less  than  the  radiation 
resistance.  Under  these  assumptions,  the  average  face  velocity  is 


avg 


<t>Kn 


R 


rad 


where  V,„  is  the  voltage  applied  to  the  electrical  terminals  and  Rrad  is  the  radiation  resistance. 
The  volume  velocity  can  be  computed  from  the  effective  area  and  the  average  face  velocity. 
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Then,  using  the  simple-source  radiation  expression,  the  ratio  of  pressure  at  one  meter  to  applied 
voltage  is 


P_ 

K 


pwfna2  IA 


R, 


rad 


In  writing  this  expression,  we  have  assumed  that  the  transducer  consists  of  two  flexural-disk 
elements. 
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Appendix  E.  Treatment  of  the  Trilaminar  Structure 

Analysis  of  a  trilaminar  flexural-disk  transducer  (see  Fig.  El)  is  a  straightforward  extension  of 
the  two-layer  transducer  with  arbitrary  neutral  plane.  In  the  trilaminar  transducer,  the  two 
ceramic  layers  are  assumed  to  have  identical  properties  and  dimensions  so  that  the  neutral  plane 
is  midway  through  the  substrate.  If  we  take  the  origin  of  the  z-coordinate  system  at  the  neutral 
plane,  then  integrations  through  thickness  of  the  ceramic  go  from  z  =  -  tc  -  tJ2  to  z  =  -  tJ2  and 
integrations  through  the  thickness  of  the  substrate  go  from  z  =  -  tJ2  to  z  =  0.  The  results  are 
then  doubled  to  account  for  the  symmetry  of  the  structure  about  the  neutral  plane. 


Fig.  El.  Trilaminar  center-supported  flexural  disk  structure.  The  two  ceramic  layers  are  assumed 
identical.  The  neutral  plane  is  in  the  middle  of  the  substrate. 


If  we  define  the  substrate  thickness  variable,  ts  to  be  half  the  actual  substrate  thickness  and  let  z0 
be  equal  to  -  ts,  then  we  can  use  the  integration-limit  constants  derived  for  the  two-layer, 
arbitrary  neutral-plane  formulation.  We  need  only  keep  in  mind  that  we  are  analyzing  only  half 
of  the  transducer  structure. 

In  equivalent-circuit  modules,  the  following  changes  are  made.  The  equivalent  mass  and 
stiffness  are  both  doubled.  The  coupling  factor  remains  the  same.  The  changes  in  other  factors 
depend  on  the  electrical  connection  of  the  two  ceramic  layers.  If  the  ceramic  layers  are 
connected  (electrically)  in  parallel,  then  the  short-circuit  charge,  the  blocked  capacitance,  and  the 
transduction  factor  are  all  doubled,  while  the  open-circuit  voltage  remains  the  same.  If  the 
ceramic  layers  are  connected  in  series,  then  the  short-circuit  charge  and  the  transduction  factor 
remain  the  same  while  the  open-circuit  voltage  doubles  and  the  blocked  capacitance  is  reduce  by 
a  factor  of  two.  For  this  release  version,  the  parallel  electrical  connection  is  assumed  (but 
conversion  of  the  parameters  for  series  connection  is  trivial). 
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Appendix  F.  Stress  Calculation  and  Interpretation 

An  important  part  of  the  design  process  is  determination  of  the  maximum  source  level  achievable 
from  the  transducer.  One  of  the  critical  limits  is  stress  in  the  ceramic.  If  the  desired  source  level 
is  prescribed,  then,  the  required  volume  velocity,  Q,  can  be  determined: 


where  p  is  the  desired  acoustic  pressure  at  one  meter,  p  is  the  density  of  the  fluid,  and /is  the 
frequency. 

Using  the  dimensions  of  the  transducer  and  the  volume  factor  from  the  shape  integral,  7V,  the 
peak  deflection  that  would  produce  that  volume  velocity  is, 

^ peak  7  7  2  r 

In  f  na  /„ 

Then,  having  the  peak  deflection  and  the  deflection  shape,  the  stresses  can  be  calculated.  In  this 
second  release,  only  the  normalized  stresses  are  calculated.  The  actual  stresses  can  be  calculated 
manually  as  outlined  above. 

The  radial  and  tangential  strains  are  determined  by  the  deflection  function,  w, 

_  d2w  _  z  dw 

\  ~  z  ,  2  >  - 

dr  r  dr 

The  stresses  can  be  determined  from  the  following  set  of  equations: 

7]  =  cc^  S ,  +  cc^2  S2  -  hh^  D3 

T2  =  ct&  S ,  +  ccfx  S2  -  hhixDi 

E3  =  -AAj,  (£,  +  S2)  +  PP33  D3 

From  the  orientation  of  the  electrodes,  Z/  is  constant  in  the  3-direction  so  we  can  integrate  the 
third  equation  over  the  thickness  of  the  ceramic  to  find  D3  in  terms  of  the  strain  and  the  voltage 
on  the  electrodes. 

The  electrode  voltage  is  calculated  by  integrating  the  electric  field  from  one  electrode  to  the 
other.  In  most  published  analyses  of  piezoelectric  transducers,  the  sign  is  handled  incorrectly. 
The  voltage  (electrical  potential)  is  minus  the  integral  of  electric  field.  (The  electric  field  is 
minus  the  gradient  of  the  potential.)  For  the  purposes  of  analysis  of  transducer  performance,  this 
error  is  generally  inconsequential;  it  can  be  “rectified”  simply  by  reversing  the  polarity  of  the 
voltage  terminals  in  the  equivalent  circuit.  However,  when  analyzing  stress  and  the  coupling 
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between  mechanical  and  electrical  properties,  the  sign  is  important.  The  first  report  did  not 
address  this  sign  error.  To  avoid  confusion  in  this  second  report,  the  terminal  voltage  will  be 
defined  so  that  a  positive  voltage  applied  to  the  electrodes  produces  an  electric  field  in  the 
positive  3-direction.  With  the  coordinate  definitions  used  so  far,  this  means  that  the  “ground” 
terminal  is  the  electrode  between  the  substrate  and  the  ceramic  and  the  “positive”  terminal  is  the 
other  electrode.  If  the  terminal  voltage  is  Vh  then 

-K  =  -\E,dz 

where  the  integration  is  taken  from  the  lower  electrode  to  the  upper  electrode.  Therefore, 

-K  =  **,.  f  (S,  *  s,)dz  -  Ws„D,t, 

-<r-*0 


Since 


5,  +  S2 


-  z 


w 


the  integration  is  straightforward: 


V,  =  ~hh 3, 


f  .  mA 

w  +  — 
r 


(t2c  +  2tc  z0 )  -  fife  Z)3  tc 


or 


D  =  Mil  (jc  +  2zo) 


W  +  — 

r 


{  rj  fi&t. 

The  maximum  stress  in  the  ceramic  is  at  the  outer  (lower)  face  where  z  =  -tc-zo'. 


ym ax 


op  max 
ll 


(fc  +  Zo) 

{tc  +  Zo)^CC12  W”  +  JCCU  W> 


f  1  \ 

cc^  w"  +  —  cc?2  w 

V  r  J 

J 


hK  A 
hK  A 


Inclusion  of  the  Di  term  in  the  stress  calculation  will  be  deferred  to  the  next  release  so  the 
stresses  calculated  in  the  second  release  should  be  considered  cautiously. 

As  an  interim  solution,  the  limiting  performance  can  be  developed  in  terms  of  maximum  normal 
strain.  This  is  not  the  preferred  criterion  for  failure  of  brittle  material  (maximum  tensile  stress  is 
the  preferred  limit)  but  it  will  serve  as  a  rough  approximation  in  this  second  release.  For  cases  of 
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practical  interest,  the  maximum  normal  strain  is  equal  to  the  radial  strain  at  the  inner  edge  (and 
outer  surface)  of  the  disk: 


sr  =  $(*  =  ->,-  z„)  = 

dr 

It  is  simpler  to  normalize  the  deflection  function,  w,  in  the  following  manner, 

^  =  w/w0  ■  1)  S  rla 

where  w0  is  the  maximum  value  of  the  deflection  (at  the  edge  of  the  center-supported  disk  or  at 
the  center  of  the  edge-supported  disk)  and  a  is  the  disk  radius.  The  maximum  strain  can  then  be 
written 


sr  =  (',  +  *,)^r(i) 

where  the  primes  denote  differentiation  with  respect  to  tj.  If  we  specify  a  maximum  strain  value, 
this  expression  permits  us  to  find  the  corresponding  peak  deflection: 


gmax 


(*c  +  *)«'(*) 


(If  we  ignore  the  stress  component  related  to  £>3,  we  could  find  the  peak  deflection  from  a 
maximum  stress  criterion.  The  first  derivative  of  £  is  zero  at  r\  =  b,  so 


wn 


a2  T™ 


<(tc  +  z0)^(b) 


for  maximum  normal  stress.) 

Having  the  peak  deflection,  the  radiated  pressure  can  be  found  from  the  simple-source 
approximation: 


v"™1  =  (O  Vf“ 

avg  ^  rvavg 


co^rT- 


where  the  shape  integrals,  /„  and  IA,  are  defined  in  the  next  section; 


=  v1™  A 

ravg 


"a1  h  C 
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and 


p(lm)  = 


for  the  pressure  at  one  meter  (assuming  the  transducer  contains  two  flexural  disk  elements). 
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Appendix  G:  Relationship  of  Solutions  to  Reference  1 

Some  adjustment  is  necessary  to  compare  results  with  Woollett01.  In  that  reference,  the 
deflection  functions  for  the  edge-supported  cases  are  normalized  so  that  the  maximum  deflection 
is  one;  however,  the  deflection  function  for  the  center-supported  case  is  not. 

The  normalized  deflection  function,  £,  used  here  always  has  a  maximum  value  of  one.  The 
deflection  function  amplitude  used  in  Ref.  1  is  a„  where  a,  is  one  for  the  edge-supported  cases 
and  a,  is  0.369  for  the  center-supported  case. 

Several  of  the  constants  used  in  Ref.  1  are  similar  to  the  shape  integrals  used  here: 


A"  =  fL 


1  P  0 


+ 


CC, 


ccx 


11/ 

D  1  P  i 


Au  = 


1  p  o 


cc 


cc} 


11/ 

D  XP\ 


_  3  gg33  hh3i  / 

A  D  \2P0 

4  cc„ 


+  h\) 


K  =  2  a]  I xz 


Using  these  associations,  the  results  from  Ref.  1  can  be  compared  to  the  present  results  by 
replacing  the  passive  substrate  used  here  with  another  ceramic  layer  having  the  same  properties 
and  dimensions  as  the  first  one.  The  comparisons  are  not  exact  because  the  basic  deflection 
functions  are  different. 


01  R.  S.  Woollett,  “Theory  of  the  piezoelectric  flexural  disk  transducer  with  applications  to  underwater  sound,”  USL 
Research  Report  No.  490,  U.S.  Navy  Underwater  Sound  Laboratory,  December  5,  1960. 
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Appendix  H:  MathCad  Code  Listing 

The  following  modules  are  listed  in  this  appendix.  Please  read  the  current  Flexlntro  file  supplied 
with  the  MathCad  files  for  updated  formats  and  program  modules. 

Flexlntro3_0.mcd 

CC2L_3_0config.mcd 
CC2L_3_0deflect.mcd 
CC2L_3_0equivckt.  mcd 

CC3L_3_Oconfig.mcd 

CC3L_3_0deflect.mcd 

CC3L_3_0equivckt.mcd 

ES2L_3_0config.mcd 

ES2L_3_0deflect.mcd 

ES2L_3_0equivckt.mcd 

ES3L_3_0config.mcd 

ES3L_3_0deflect.mcd 

ES3L_3_0equivckt.mcd 


Performance3_0.mcd 
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Flexural-Disk  Transducer  Analytical  Model  -  DRAFT  VERSION  3.0 
Introduction 

This  model  presents  an  analytical  solution  for  the  flexural-disk  transducer  in  a  manner  similar  to  that 
used  by  Woollett  (Theory  of  the  Piezoelectric  Flexural  Disk  Transducer  with  Applications  to  Underwater 
Sound,  USL  Report  490,  December  5,  1960).  The  primary  focus  is  on  transmitting  transducers  with 
center-supported  disks  but  edge-supported  disks  and  receiving  performance  will  also  be  treated. 

The  three  most  apparent  differences  between  this  work  and  Woollett's  report  are:  (1)  higher-order 
deflection  functions  are  used  to  describe  operation  in  the  vicinity  of  resonance,  (2)  a  non-zero  diameter 
center  support  is  considered,  and  (3)  arbitrary  location  of  the  neutral  plane  is  permitted.  Allowing  a 
non-zero  center  support  diameter  is  critical  in  characterizing  the  stress  distribution.  Arbitrary  location  of 
the  neutral  plane  permits  modeling  of  two-layer  (ceramic/substrate)  disks. 

Users  should  be  aware  that  this  is  a  draft  version  rather  than  a  fully  validated  production  release. 


Transducer  Types 

Four  configurations  of  the  flexural-disk  transducer  are  analyzed.  Two  configurations  are 
center-supported.  The  first  center-supported  configuration  consists  of  a  single  ceramic  layer  and  a 
substrate  layer;  the  second  consists  of  a  substrate  layer  and  two  ceramic  layers  on  either  side  of  the 
substrate.  These  configurations  are  shown  below: 
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Two  configurations  are  edge-supported.  The  first  edge-supported  configuration  consists  of  a  single  layer  of 
ceramic  and  a  substrate  layer;  the  second  consists  of  a  substrate  layer  and  two  ceramic  layers  on  either 
side  of  the  substrate.  In  each  case,  the  composite  disk  is  simply  supported  at  the  edge: 


HINGED  EDGE 
SUPPORT 


HINGED  EDGE 
SUPPORT 


Third  Release 


In  this  third  release,  the  elements  in  the  two-layer  structures  are  considered  to  be  two-layer  laminations 
of  ceramic  and  substrate  with  arbitrary  properties  and  thicknesses.  Consequently,  the  neutral  plane  is  not 
restricted  to  being  located  at  the  ceramic-substrate  interface.  Some  caution  must  be  exercised  now  in  that 
it  is  possible  to  put  the  neutral  plane  in  the  ceramic,  which  can  cause  failure  under  hydrostatic  load  if  the 
structure  is  not  prestressed.  Also,  in  this  second  release,  an  approximate  method  for  calculating  the  mode 
shape  and  shape  integrals  is  used.  Incorporation  of  this  energy-based  approximate  method  will  eventually 
permit  modeling  of  structures  in  which  the  ceramic  does  not  completely  cover  the  substrate.  Version  3.0 
calculates  the  equivalent  mechanical  resistance  and  the  dielectric  loss  as  equivalent-circuit  elements. 

There  are  three  important  additions  in  version  3.0.  First,  the  three-layer  edge-supported  disk  structure  is 
treated  (ES3Lxxx).  Second,  a  feature  has  been  added  to  the  Performance  Module  that  permits 
incorporation  of  a  fluid-filled  cavity  behind  the  flexural  disk.  Finally,  physical  units  have  been  incorporated  in 
the  Performance  module.  Since  files  do  not  preserve  the  units  for  quantities,  the  calculations  in 
Performance  were  done  previously  without  units  (but  assuming  SI  units).  In  this  version,  the  units  are  now 
explicit  (and  SI). 

As  in  the  previous  version,  this  release  is  presented  as  a  series  of  modules.  This  architecture  is 
described  below  and  is  convenient  for  exploring  parameter  variations. 


Assumptions 

A  number  of  assumptions  have  been  made  in  this  version  of  the  flexural  disk  model.  Several  of  these 
assumptions  will  be  relaxed  as  the  model  is  developed;  some  may  be  retained  in  order  to  preserve  the 
analytical  nature  of  the  model. 

Thin-plate  theory  (plane  stress;  no  shear) 

Ceramic  layer  extends  over  entire  diameter  of  disk  structure 

Simple  assumed  deflection  functions  are  not  used  as  there  is  no  significant  saving  in  computation  nor  is 
there  any  reduction  in  the  analytical  nature  of  the  model  by  using  "exact"  mode  functions  or  mode  functions 
developed  by  Rayleigh-Ritz  from  higher  order  polynomials.  When  calculation  of  hydrostatic  stresses  are 
added  to  the  model,  the  exact  static  deflection  will  be  used  for  the  same  reasons. 


Organization 

The  naming  convention  for  the  modules  is  as  follows.  The  first  part  of  the  name  indicates  the  transducer 
configuration.  CC2L  is  the  center-supported,  two-layer  configuration;  CC3L  is  the  center-supported, 
three-layer  configuration;  ES2L  is  the  edge-supported,  two-layer  configuration;  and  ES3L  is  the 
edge-supported,  three-layer  configuration.  For  each  transducer  type,  there  are  three  modules  and  these  are 
indicated  by  the  second  part  of  the  module  name.  The  first  module  is  the  configuration  module 
(xxxx_config),  the  second  (xxxx_deflect)  is  the  module  that  finds  the  mode  deflection  function  and  the 
shape  integrals,  while  the  third  module  (xxxx_equivckt)  computes  the  equivalent-circuit  parameters  and 
some  of  the  important  operating  parameters. 

Run  one  of  the  configuration  modules  (CC2L_config,  CC3L_config,  ES2L_config,  or  ES3L_config)  first  to 
define  the  transducer  configuration,  the  piezoelectric  material,  and  the  substrate  material.  Many  of  the 
parameters  subsequently  calculated  do  not  depend  on  specific  dimensions.  In  those  cases,  dummy 
dimensions  can  be  entered.  The  first  release  retained  sufficient  flexibility  to  calculate  impractical  limiting 
cases  (such  as  the  center-hole  radius  going  to  zero).  Unrealistically  small  center  radii  should  not  be  used 
in  this  version. 


The  configuration  module  also  reads  the  piezoelectric  parameters  and  substrate  properties  from  the 
materials  files  and  derives  the  necessary  plane-stress  forms  of  other  parameters  for  subsequent 
calculation.  Manufacturers'  tabulated  values  MUST  NOT  BE  USED  for  the  c  or  h  parameters  or  for  the 

clamped  dielectric  permittivity,  0s.  Tabulated  values  are  not  appropriate  for  2D  plane-stress  problems. 

Then  run  one  of  the  deflection-function  modules  (CC2L_deflect,  CC3L_deflect,  ES2L_deflect,  or 
ES3L_deflect).  This  module  finds  the  first  mode  eigenvalue  and  the  corresponding  deflection  shape.  Then 
the  so-called  shape  integrals  -  factors  that  depend  only  on  the  shape  of  the  deflection  curve  --  are 
computed.  The  shape  integrals  are  used  in  the  calculation  of  equivalent-circuit  parameters  and  other 
performance  parameters. 

Finally,  run  the  parameter  module  (CC2L_equivckt,  CC3L_equivckt,  ES2L_equivckt,  or  ES3L_equivckt) 
to  evaluate  the  equivalent-circuit  parameters  for  the  transducer  configuration  specified  in  the  configuration 
module.  Once  the  equivalent-circuit  parameters  have  been  stored,  Performance  can  be  run  to  compute 
admittance,  TVR,  TCR,  and  a  few  other  operating  properties.  The  inputs  to  set  up  the  optional  fluid-filled 
cavity  are  entered  in  Performance. 


File  Structure 

Piezoelectric  materials  files  (PZT4.prn,  PZT5H.prn,  PZT8.prn): 

Each  piezoelectric  material  file  contains  seven  values  from  various  manufacturer  data  sheets.  The  seven 
values  are: 

s°i2  s°i2  ^T33^o  d3i  □MIlanDD  Qmech 

Files  for  other  materials  can  be  generated  easily.  Create  an  ASCII  file  (use  the  extension  .prn)  with  these 
seven  values  entered  on  a  single  line  separated  by  spaces.  The  values  are  stored  in  the  file  as  they  are 
usually  listed  in  properties  tables:  the  compliances,  s,  in  units  of  10-12  m2/N;  the  dielectric  permittivity  as 
relative  dielectric  permittivity;  and  the  d  coefficient  in  units  of  1CM2  mA/.  These  values  are  the  same  for  3D 
and  2D  problems  so  they  may  be  transcribed  directly  from  manufacturer  data  sheets.  The  configuration 
module  reads  these  values  and  then  multiplies  the  compliances  and  the  d  coefficient  by  1CH2  and 
multiplies  the  relative  dielectric  permittivity  by  the  permittivity  of  free  space. 

Substrate  properties  files  (brass. prn,  aluminum.prn,  steel.prn): 

The  substrate  properties  files  each  contain  three  values:  the  density  in  kg/m3,  the  elastic  modulus  in  GPa, 
and  the  Poisson's  ratio.  The  user  can  generate  other  substrate  files  as  ASCII  files  (with  .prn  extension)  as 
described  above. 

Configuration  files  (xxxx_config.prn): 

The  configuration  files  contain  15  values  -  three  dimensions  and  seven  physical  properties  of  the  ceramic, 
three  properties  and  the  thickness  of  the  substrate,  and  the  location  of  the  neutral  plane: 


□c  a  b  tC  °D11  c°12  h31  ^33  tan0  Qmech  Psub  Esub  vsub  ‘sub  z0 


Units  are  not  retained  in  the  write  operation  so  all  values  are  converted  to  proper  SI  units  before  being 
written.  User  input  values  can  be  entered  in  other  systems  as  long  as  they  are  entered  with  those  units.  If 
no  units  are  specified,  SI  units  will  be  assumed.  The  outer  radius  of  the  disk  is  a;  the  thickness  is  tc;  and 
the  ratio  of  inner  radius  to  outer  radius  is  b.  For  the  edge-supported  disk,  b  is  set  to  zero  in  the 
configuration  module.  The  density  and  the  loss  tangent  are  transferred  directly  from  the  properties  file.  The 
stiffnesses,  cD,  the  h  coefficient,  and  the  clamped  permittivity  are  derived  based  on  the  plane-stress 
assumption  appropriate  to  the  plate  bending  equations.  These  values  will  not  compare  to  published  values. 


Shape  files  (xxxx_deflect.prn): 

The  shape  files  hold  the  values  for  the  shape  integrals;  the  peak  deflection,  wnmax,  for  unity  strain;  the 
size,  N,  of  the  coefficient  matrices;  and  the  polynomial  coefficient  matrix,  CC: 


'po 


■pi 


'ke 


wnmax  N  CC 


The  shape  integrals  are  functions  only  of  the  shape  of  the  deflection  function  and  so  permit  rather  general 
conclusions  to  be  drawn  about  different  configurations  without  substituting  specific  dimensions  or  physical 
properties.  Their  value  is  not  fully  exploited  by  this  initial  draft.  The  first  two,  lp0  and  Ip1,  are  related  to 
elastic  strain  energy.  The  third,  iq,  is  related  to  the  charge  produced  under  strain.  The  fourth,  lv,  relates  the 
average  face  velocity  to  the  peak  face  velocity.  The  fifth,  lKE,  is  related  to  the  kinetic  energy.  And  the 
sixth,  ia,  is  related  to  the  active  face  area.  The  shape  integrals  and  the  configuration  files  are  used  to 
compute  the  equivalent-circuit  parameters  and  various  performance  measures.  The  peak  deflection  for  unity 
strain  can  be  used  to  calculate  performance  limits  (this  is  implemented  in  the  third  module  for  each 
configuration).  N  and  CC  can  be  used  to  construct  the  actual  deflection  function. 


Circuit  parameter  files  (xxxx^equivckt.prn): 

The  equivalent-circuit  parameter  files  hold  the  following  values: 


mass  stiffness  0  CEB  k2  tan<5  Rmech  r_rad  m_rad  Aeff 

These  parameters  are  the  equivalent  mass  and  stiffness  at  resonance;  the  electromechanical  transduction 
factor,  0;  the  blocked  electrical  capacitance,  CEB;  the  coupling  factor,  02;  the  loss  tangent,  tanO,  of  the 
ceramic;  and  the  equivalent  mechanical  resistance,  Rmech,  of  the  composite.  The  quantites,  r_rad  and 
M„rad  are  the  radiation  resistance  and  radiation  mass  for  the  parallel  equivalent  for  radiation  impedance. 
Aeff  is  the  effective  piston  area. 

Results 


Mode  deflection  function 

Basic  energy  quantities 

Unloaded  and  loaded  resonance  frequency 

Coupling  factor 

Equivalent-circuit  parameters 

Admittance  (unloaded  and  loaded) 

Transmitting  voltage  response  (TVR) 

Transmitting  current  response  (TCR) 

Low-frequency  free-field  voltage  sensitivity  (FFVS) 
Calculation  of  maximum  stress  in  ceramic  and  substrate 
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Flexural-Disk  Transducer  Analytical  Model  -  DRAFT  VERSION  3.0 


This  worksheet  is  specific  to  the  annular  (center-supported)  disk.  First,  the  physical  dimensions  are 
specified  and  the  piezoelectric  properties  file  is  read.  Then,  the  two-dimensional  properties  are  derived. 
Finally,  the  results  are  written  to  a  configuration  file. 


I.  SELECT  MATERIALS  [Highlighted  values  are  user  inputs] 

Choose  from  PZT4,  PZT5H,  or  PZT8  for  file  name  in  pzt  READPRN. 
pzt :  =  READPRN(  PZT4 ) 


sell  :=pztM-10T12-J- 


se!2  :=pzt  .*10" 


et33  =  pztQ  2-e0 


d31  :=pzt  •  10'12— ! — 
°*3  volt 


P  :  =  pzt 


0,4  3 

m 


tan5  :  =  pzt 


Qmech  =  pzt 


Choose  from  steel,  aluminum,  or  brass  for  file  name  in  substrate  READPRN. 
substrate  :=  READPRN( brass) 


9 

Esub  :=substrateQ  ^10  -Pa  vsub  :=substrateQ  2 

Esub  =  1.04-1 011  -kg-m-|-sec"2 


vsub  =0.37 


//.  SELECT  DIMENSIONS 


Enter  dimensions  of  ceramic  disk  (a  =  outside  radius;  bi  =  inner  radius;  tc  =  thickness) 
and  thickness  (tsub)  of  substrate  disk: 

[Enter  units  for  in,  mm,  cm,  or  m;  if  no  units  are  entered,  meters  will  be  assumed] 

a  =7-cm  bi:=l-cm  tc:=2.5-mm  tsub:  =  4-mm 

a=0.07-m  b  =0.143  tc=2.5-l(f3  -m 


III.  DERIVE  PROPERTIES  FOR  TWO-DIMENSIONAL  ANALYSIS 


sdl  1  :=sel  1  - 


d3 1 2 
et33 


d3 1 2 

sdl 2  :  =  sel2 - 

et33 


sdl  2 
sdl  1 


0.485 


cdl  1 


sdl  1 

sdl 1 2  —  sdl22 


cdl2  :  = 


-  sdl2 

sdl l2  -  sdl22 


cdl  2 
cdl  1 


=  0.485 


es33  :=et33  - 


2-d3 1 2 
sel  1  f-  sel2 


h31 


d31 

st33-(sdl  1  -i-  sd!2) 


IV.  CALCULATE  LOCATION  OF  NEUTRAL  PLANE 

Origin  of  z-coordinate  at  neutral  plane;  zO  gives  distance  from  neutral  plane  to  ceramic  substrate 
interface: 


fs  -  - 


Esub 

2 

1  -  vsub 


zO 


tsub2-fs-  tc2*  cd  1 1 
2*(tsub-fs-t-  tc  cdl  1 ) 


zO  =  7.6051 10  4  -m 


WARNING:  If  zO  is  negative,  then  the  neutral  plane  is  in  ceramic  and  some  part 
of  the  ceramic  will  be  in  tension  under  hydrostatic  load  unless  prestressed. 

V.  WRITE  CONFIGURATION  FILE 

PRN  files  cannot  handle  dimensions  so  all  quantities  are  written  (in  SI)  without  units. 


out, 


0,4  ■' 


cdl  1 
cO 


out 


0,1 


out 


0,5 


m0 

cdl  2 
cO 


OUt0,2;  =  b 


out 


0,6  ' 


h31 

hO 


out, 


0,3 


tc 

m0 


out 


0,7 


ss33 

eO 


outQ  :  =  tan8  outQ  g  =Qmech 


out 


o,io  ’ 


psub 
pO 


out, 


0,11 


Esub 

cO 


out 


0,12 


:=  vsub 


out 


0,13 


tsub 

mO 


out 


0,14 


zO 

mO 


WRITEPRN(  CC2L_config)  -  out 


The  configuration  file  contains  the  following  quantities  in  this  order: 

density  of  ceramic,  p 

outer  radius  of  disk,  a 

ratio  of  inner  radius  to  outer  radius,  b 

thickness  of  ceramic,  t 

cdll  for  ceramic 

cd12  for  ceramic 

h31  for  ceramic 

es33  for  ceramic 

loss  tangent  for  ceramic,  tan5 

mechanical  Q  of  ceramic 

density  of  substrate,  psub 

modulus  of  substrate,  Esub 

poisson's  ratio  of  substrate,  vsub 

thickness  of  substrate,  tsub 

z  location  of  neutral  plane,  zO,  with  respect  to  bottom  of  ceramic 


Unit  normalization  quantities  (predefined): 

p0  =  l-—  m0=l-m  cO=l-Pa 

3 
m 


h0=  1- 


volt 


m 


n-i  farad 

e0  =  1 - 

m 


£0  =  8.854-10 


.  12  farad 


m 
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Flexural-Disk  Transducer  Analytical  Model  ■  DRAFT  VERSION  3.0 

Module  2.CC2L  deflect:  Center-Supported  Disk  -  Clamped  Inside 

This  worksheet  reads  the  configuration  file,  calculates  the  first  axisymmetric  mode  deflection  function, 
and  then  calculates  all  of  the  shape  integrals.  The  shape  integrals  are  written  to  a  file  for  subsequent  use, 
In  this  version,  the  polynomial  Rayleigh-Ritz  method  is  used  to  compute  the  mode  deflection  function.  An 
eighth-order  polynomial  is  used,  which  produces  accurate  results  in  all  cases  except  unrealistically  small 
inner  radii.  This  module  is  expecting  the  configuration  file  from  version  2  -  that  is,  the  configuration  with 
arbitrary  location  of  the  neutral  plane  and  specific  properties  for  the  substrate. 


I.  READ  CONFIGURATION  FILE 

(pc  a  b  tc  edit  cd!2  h31  es33  tanS  Qmech  psub  Esub  vsub  tsub  zO )  :=  READPRN(CC2L_config) 

a  -0.07  b  =0.1429  tc  =0.0025  vc:  =  -^^  vc  =0.485282 

cdl  1 

v  :=  0.5-(vc  +-  vsub)  v  =  0.427641 

II.  MODE  SOLUTION 

In  this  section,  the  Rayleigh-Ritz  technique  is  used  to  find  the  fundamental  eigenvalue  and  eigenfunction  for 
an  8th  order  polynomial.  (N  is  preset  to  8  below.) 

UU  :=U(N,b,v)  TT  :  =  T(N,b,v) 

gvats  :  =  genvals(UU,TT) 

gvecs  genvecs(UU,TT) 

sgcomb  :  =  rsort (stack (gvalsT,  gvecs)  ,o) 

fundval  :  =  Jsgcomb0  Q  fundval  =4.692896 


III .  FIND  THE  COEFFICIENTS  FOR  THE  DEFLECTION  FUNCTION 
fundcoef :  =  submatrix(  sgcomb,  1 ,  N ,  0 , 0 ) 

CC  :  =  scaledcoef(  fundcoef,  N,  b ) 


The  coefficients  are  normalized  so  that  the  maximum  deflection  is  one. 


IV.  PLOT  THE  NORMALIZED  STRESSES 


eta  :=b,b  +■  0.002..  1 


V.  CALCULATE  THE  SHAPE  INTEGRALS 


Ike  :=CCT  TT-CC 

Ike  =  0.206775 

la  :  =  1  -  b2 

la  =0.97958 

Iv  :  =  CCT’VV(N,b) 

Iv=  0.569635 

Iq:=CCTQQ(N,b) 

Iq  =  1.221375 

Ipsum  :  =  CCT-Usum(N,b,v)-CC 

Ipsum  =  5.407692 

Ipl  :=CCT-Upl(N,b, v)-CC 

Ipl  =  1.491756 

IpO  :=  Ipsum  -  Ipl 

IpO  =3.915936 

The  maximum  strain  is  at  the  inner  radius.  The  second  derivative  of  the  normalized  deflection  is 
calculated  below  (maxcurve)  and  the  peak  deflection  (i.e.,  deflection  at  r  =  a)  corresponding  to  a  strain  of 
one  is  saved  (wnmax). 

2 

maxcurve  :  =  CCT-SM(N,b)  wnmax  :  = -  wnmax  =0.136099 

(tc-r  zO)  maxcurve0 


Write  shape-integral  file: 

WRITEPRN(CC2L_deflect)  :=(lp0Q  IplQ  IqQ  IvQ  IkeQ  la  wnmax  N  CC) 

Because  several  of  the  shape  integrals  are  calculated  through  matrix  operations,  they  remain  as 
matrices  even  though  having  dimension  1x1 .  The  subscript  index  is  used  to  convert  to  scalar  form  before 
writing.  The  deflection  coefficients  are  written  as  a  matrix  (Nxl)  as  CC. 


What  follows  are  predefinitions... 


N=8 


scaledcoef(  vect ,  N,  b )  = 


peakdefl  <-  0 
for  iie  1  ..N 


•  4  /  *  •  4  V  1  11  •  •  t  11  I 

incr<-  1  -  (ih-  1  )b  +-  n-b 
peakdefl  peakdefl  +  incr*  vect..  ] 


vectscaled<— 


vect 

peakdefl 


vectscaled 


Usum(N,b,v)  = 


for  me  1 .. N 
for  n  e  1 ..  N 

i  _ 

ul<—  (m-h  1  )-(n  +■  1) - 

m  i-  n 


u2*~- (in-  1  )-b - +  -(mr  l)-b - 

n  m 

uam_  i  n_  1  )(n  +-  1  )-(ul  +  u2-h  -  bm_hn-ln(b)) 


ua 


Upl(N,b,v)s 


for  me  1 ..  N 
for  ne  1 ..  N 


ubb  .  n*(n  +  l)*(m+-  !)• 

m  -  1 ,  n  -  1  x  '  v  7 


l-bm+n  bm-(l-bn) 


nun 


ub<~  ubb  -h  ubb 
ub 


U(N,b,v)s  uu<—  Usum( N, b ,  v)  -  ( 1  -  v)-Upl(N,b,v) 
uu 


SM(N,b)  = 


for  me  1  ..N 
sm  <-m-(m-h  1  )-bm~  1 

m-  1  v  7 


sm 


T(N,b,  v)  = 


for  me  1  ..N 
for  ne  1  ..N 


tU 


1  -  b 


m  -f  n-f  4 


m  +  n  -1-  4 


(n-t-  l)  b 


n  1  —  b 


m  -(-4 


m-f-  4 


t2<-nb 


n-f- ! 


1  -  b 


m-f  3 


(m-t-  l)-b 


1  -  b' 


n-f  4 


m-f  3  n  -f  4 

[4 


t3<—  (m-f  l)-(nf  l)-bm  — —  -  (2-m-ntmtn)'bmi‘nf1'- 

4 

^  Lmf  1  1  -  bn^  m-fn-f2  1  -  b 

t4<—  m  b  ^ - -fm-n-b  ^  ^ - 

n  f  3  2 

tt  .  <— tl +- 12 -f  t3 -r  t4 

m -  1  ,n-  l 


VV(N,b)  = 


for  me  1  ..N 


matv 


m~  1 


1  -  b 


m-f  3 


m-f  3 


-  (m-f  1  )  b 


l-bJ 


-f  m  b 


m  -f  1  1  “  b 


2*  matv 


QQ(N,b)- 


for  me  1  ..N 

matvm_  (m -f  1  )*(l  -  bm) 


matv 


slope(  eta,  vect,  N,b)  = 


slopes- 0 
for  iie  1 ..  N 

slope*- slope  -f  vect.._  ^(ii  -f  1  )-(eta11  -  b11) 


slope 


curve(  eta ,  vect ,  N ,  b )  = 


curve*—  0 
for  iie  1  ..N 

curve*-  curve  -f  vect.._  *(ii  -f  1  J-ii-eta1 


ii-  l 


RadialStress( eta,  vect ,  N, b ,  v)  =  curve( eta,  vect ,  N,  b )  -f 


v*  slope(  eta,  vect ,  N,  b ) 


eta 


TangentiaIStress(  eta,  vect ,  N,  b ,  v)  =  v-  curve(  eta,  vect ,  N,  b)  ± 


slope(eta,vect,N,b) 

eta 


F£NNSt*T£ 
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Flexural-Disk  Transducer  Analytical  Model  -  DRAFT  VERSION  3,0 

Module  3:  Equivalent-Circuit  Parameters  -  Bilaminar  Flex  Disks 

This  worksheet  reads  the  configuration  and  shape-integral  files  for  a  particular  case  and  computes  the 
equivalent-circuit  parameters. 


V 

avg 


In  these  equivalent  circuits,  the  mechanical  quantities,  force  (pA  or  acoustic  pressure  times  effective  area) 
and  velocity,  v,  are  related  to  the  electrical  quantities,  voltage,  e,  and  current,  i.  CEB  is  the  blocked 
electrical  capacitance  and  RD  is  the  dielectric  loss  (tan5/w*CEB).  The  mechanical  mass,  stiffness,  and 
mechanical  damping  are  m,  k,  and  Rm  respectively.  The  mechanical  resistance  is  computed  based  on  the 
mechanical  Q  of  the  ceramic  and  the  fraction  of  strain  energy  that  is  stored  in  the  ceramic  layer  of  the 
composite.  The  radiation  impedance,  Zrad,  is  represented  as  a  parallel  combination  of  resistance  and 
mass  (see  below).  Finally,  the  electromechanical  transduction  factor  (the  "turns  ratio")  is  q. 

I.  READ  FILES 
Read  configuration  file: 

(pc  a  b  tc  cdll  cdl2  h31  es33  tanS  Qmech  psub  Esub  vsub  tsub  zO )  :=  READPRN(CC2L_config) 
a  =0.07  b  =0.143  tc=2.5-10~3 


Read  shape-integral  file: 


(IpO  Ipl  Iq  Iv  Ike  la  wnmax  N  CC)  :-READPRN(CC2L_deflect) 

Aeff  :  =  7i*a2*Ia  Aeff  =  0.015 

zO 
tsub 


acl  •—  1  "h  2- 


zO 

tc 


zO 


fzO 


ac2  :=  1  -h  3 —  +  3 

tc  \tc 


as2  : 


1  -  3—^-  -t-  3- 
tsub 


II.  CALCULATE  EQUIVALENT-CIRCUIT  PARAMETERS 
Equivalent  mass: 

la2 

mass  :=2-7i‘a2-Ike*— -(pc-tc-h  psub  tsub)  mass  =0.998 

Iv2 


Equivalent  stiffness: 

.2 

fll  :=cdll  -  0.75-es33-h312--— - 

ac2 


fl2  :~cdl2  -  0.75  es33-h312  — — 

ac2 


csl  1 


Esub 
1  -  vsub2 


cs!2  :=  vsub-csl  1 


Iff fl  MpO+  fl 2*Ip  1 


Iss  =  csl MpO-t-  cs  1 2- Ip  1 


stiffness  :  =  — -(tc3-ac2-Iff-f-tsub3-as2-Iss)  stiffness  =4.157*  IO7 

_  2  T  2 
3*a  Iv 


Transduction  factor: 


q  :  =  -  7i-h31-es33Tq- 


— -tc-acl] 

iv  ; 


<|)  =0.392 


Blocked  electrical  capacitance: 

2 

CEB  :=  £s33'71' -  -la  CEB  =4.766- 10' 


tc 


Coupling  factor: 


Iv 

kfl  :=  stiffness - 

la2 


k2  := 


kf2 

kfl  -h  kf2 


kf2  :  =  jc-  es33-  h3 1 2-tc3  • 

la 


k2  =0.072 


Calculate  equivalent  Q  for  composite  structure: 


Uratio  :  = 


(tsub\3  as2  Iss 
tc  /  ac2  Iff 


Qeff  :-Qmech-(  1  -t  Uratio) 


Rmech 


stiffness  mass 
^  mass  Qeff 


Calculate  circuit  elements  for  parallel  form  of  radiation  impedance.  This  form  is  preferred  over  the 
series  form  since,  in  the  parallel  form,  the  resistance  element  is  frequency  independent.  Furthermore,  the 
parallel  form  degrades  more  gracefully  for  ka  approaching  (or  greater  than)  one.  Also,  the  assumption  is 
made  that  these  transducer  elements  would  always  be  used  in  pairs  placed  back-to-back  and  driven  in 
phase.  In  this  drive  mode,  the  radiation  impedance  is  equivalent  to  that  of  a  single  element  in  an  infinite 
rigid  baffle. 


^rad 

I T 


p_water  :  =  1 000  c_water :  =  1 500 

m_rad  :  =— -pjwater-a- (n-a^-Ia  R_rad  :=  1.44-p_water*c_water (n-a2)*Ia 

3*71 

Write  equivalent-circuit  parameters  to  a  file: 

WRITEPRN(CC2L_equivckt)  :  =  (mass  stiffness  0  CEB  k2  tanS  Rmech  R_rad  m_rad  Aeff) 


III.  MISCELLANEOUS  CALCULATIONS 


Unloaded  resonance  frequency: 


mech_res  — —  jst^ncss_  mech_res  =  1.027*  103 

2-n  'V  mass 

Rough  estimate  of  water  loading  (piston  with  uniform  face  velocity): 


rad_mass  — p_water-7i;-a3Ta 
3-tc 


co2water  := 


stiffness 

mass  -h  rad„mass 


mech_res_water 


1 

2-n 


, i 


co2water 


mech_res_water  =  745.628 


p_water-a 

pC'2-tc 


1.842 


wavenum 


2-n-  mech_res_water 
c_water 


ka  :=  wavenuma  ka  =0.219 


radR  :=  p_water-c_water-7i;*a  -la -  radR  =  540.586 

2 

Estimated  maximum  source  level  (if  strain  limited): 

[Note:  Maximum  strain  is  not  an  appropriate  failure  criterion  for  a  ceramic  material.  It  is  used  in  this 
draft  version  as  a  crude  estimate  of  limiting  performance.  Eventually,  the  more  appropriate  maximum 
tensile  stress  criterion  will  be  implemented.  The  often-used  von  Mises  stress  criterion  is  also  inappropriate 
for  failure  in  ceramic.] 

Enter  maximum  allowable  strain  Smax  0.001 
2 

Qmax  :  =  7t-a  ^-TC-mech^res^waterwnmax-SmaxTv 


pi  :  =  p_water*mech_res_water*Qmax  SL  :  =  20-log(pl )  +  120 

Maximum  pressure  at  1  meter  (Pa):  pi  =4.169*  103 


Maximum  SL  in  dB  with  respect  to  1  micropascal  at  one  meter:  SL  =  192.4 


Estimated  peakTVR: 

What  follows  is  an  optimistic  estimate  of  the  peak  transmitting  voltage  response;  mechanical  loss 
in  the  ceramic  is  ignored. 

2  _  <b 

tvr_peak  :=  p  watermech_res_water7t-a  Ta - 

radR 

TVR  in  pascals  at  one  meter  per  volt  input:  tvr„peak  =  8. 1 59 

TVR  in  dB  with  respect  to  1  micropascal  at  one  meter  per  volt:  20-log(  |tvr_peak| )  -j-  120  =  138.2 


Low-frequency  receive  response: 


The  receiving  response  is  not  normally  important  but  can  be  useful  in  reciprocity  calculations. 

ffvs;=£i!j£^. - ! - \ 

<h  .  CEB -stiffness 

1  + - 2 - 

\  <t>  / 


FFVS  in  open-circuit  volts  per  pascal:  ffvs  =  2.771*  10 

FFVS  in  dB  with  respect  to  one  volt  per  micropascal:  20*log(  |ffvs| )  -  120  =-171.1 


The  quantity  in  square  brackets  in  the  equation  for  ffvs  should  equal  02: 


- : - =  0.072  k2  =  0.072 

.  CEB-stiffness 
1  h- - 


2 
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Module  1.CC3L  confiq:  Configuration  Definition  -  Center-Supported  Trilaminar 
Disk 

This  worksheet  is  specific  to  the  annular  (center-supported)  disk.  First,  the  physical  dimensions  are 
specified  and  the  piezoelectric  properties  file  is  read.  Then,  the  two-dimensional  properties  are  derived. 
Finally,  the  results  are  written  to  a  configuration  file. 

* — - 1 


b*aK-j 


CERAMIC 

gj 

£l 

SUBSTRATE  ^ 

rs  xy 

r— 

•-  :  r,^1 

ll 

CERAMIC  * 

< 

CENTER 

tc 

SUPPORT  POST 

I.  SELECT  MATERIALS  [Highlighted  values  are  user  inputs] 

Choose  from  PZT4,  PZT5H,  or  PZT8  for  file  name  in  pzt  READPRN. 
pzt :  -  RE ADPRN(  PZT4 ) 

sell  :=pzt0  0-10'12~ 

d31  :=pzt  -lO12-^- 
0,3  volt 

Choose  from  steel,  aluminum,  or  brass  for  file  name  in  substrate  READPRN. 

c 

substrate  :=  RE ADPRN( brass) 

psub  :=  substrate0  Q-—  Esub  :=  substrate0  109Pa  vsub  :=  substrate^  2 

m3 

psub  =  8.5,103  'kg-m-3  Esub  =  1.04'  1011  -kg^m  1 ‘sec”2 


sel2  :=pzt  ,10"12-— 
°’ 1  Pa 


et33  -=pzt0  2-s0 


P-Pzt0>4-^ 

m 


tan8  :=pzt 


0,5 


Qmech  :  =  pzt 


0,6 


vsub  =0.37 


II.  SELECT  DIMENSIONS 


Enter  dimensions  of  ceramic  disks  (a  =  outside  radius;  bi  =  inner  radius;  tc  =  thickness) 
and  thickness  (tsub)  of  substrate  disk: 

[Enter  units  for  in,  mm,  cm,  or  m;  if  no  units  are  entered,  meters  will  be  assumed] 


a  =  4.0-in 


a  =  0.102‘m 


bi  :=0.25-in 


b  =0.063 


tc  :=  0.10-in 


tsub  :=  0.10- in 


tc  =  2.54*  10  *  -m 


ts2  := 


tsub 


NOTE:  For  the  trilaminar  configuration  it  is  assumed  that  the  two  ceramic  disks  are  identical. 
III.  DERIVE  PROPERTIES  FOR  TWO-DIMENSIONAL  ANALYSIS 


sdl  1  :=  sell  - 


d3l" 


sd!2  :=se!2 - 


d3K 


cdl 1 


8133 

sdll 


sdl l2 -  sdl22 


cdl2  :  =  - 


st33 
-  sdl 2 

sdll2-  sd!22 


sdl  2 
sdll 

cdl2 
cdl  1 


-  0.485 


=  0.485 


es33  :=et33 


2*d31 


h3 1 


d31 


se  1 1  -r  se  1 2  et33«(sdl  1  -h  sdl2) 

IV.  WRITE  CONFIGURATION  FILE 

PRN  files  cannot  handle  dimensions  so  all  quantities  are  written  (in  SI)  without  units. 


out 


<4  3. 

II 

o 

o 

0UV.~ 

OUt0,2:  =  b 

OUt0,3 

tC 

mO 

cdl  1 

°'4  cO 

cdl  2 

outo  s 

0>5  cO 

.  ■  h31 
°’6  hO 

OUt0,7 

_  es33 

eO 

:o,8:=tan5 

outQ  9  :-Qmech 

psub 

JUt°’10’  p0 

OUt0,l. 

Esub 

cO 

OUt0  12  '=  VSUk 


out 


0,13 


,  ts2 
mO 


out, 


0,14 


.  ts2 
mO 


WRITEPRN(  CC3L_config )  :  =  out 


The  configuration  file  contains  the  following  quantities  in  this  order: 

density  of  ceramic,  p 

outer  radius  of  disk,  a 

ratio  of  inner  radius  to  outer  radius,  b 

thickness  of  ceramic,  t 

cdll  for  ceramic 

cd12  for  ceramic 

h31  for  ceramic 

es33  for  ceramic 

loss  tangent  for  ceramic,  tan§ 

mechanical  Q  of  ceramic 

density  of  substrate,  psub 

modulus  of  substrate,  Esub 

poisson's  ratio  of  substrate,  vsub 

half  thickness  of  substrate,  ts2 

location  of  neutral  plane  (same  as  half  substrate  thickness) 


Unit  normalization  quantities  (predefined): 

pO  =  1  .iSS-  mO  =  1  •  m  cO^l-Pa 

3 
m 


h0=  1 


volt 


m 


farad 

e0=  l - 

m 


eO  =  8.854-10' 


-12  farad 


m 
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Module  2.CC3L  deflect:  Center-Supported  Trilaminar  Disk  -  Clamped  Inside 

This  worksheet  reads  the  configuration  file,  calculates  the  first  axisymmetric  mode  deflection  function, 
and  then  calculates  all  of  the  shape  integrals.  The  shape  integrals  are  written  to  a  file  for  subsequent  use. 
In  this  version,  the  polynomial  Rayleigh-Ritz  method  is  used  to  compute  the  mode  deflection  function.  An 
eighth-order  polynomial  is  used,  which  produces  accurate  results  in  all  cases  except  unrealistically  small 
inner  radii.  This  module  is  expecting  the  configuration  file  from  version  2  -  that  is,  the  configuration  with 
arbitrary  location  of  the  neutral  plane  and  specific  properties  for  the  substrate. 


/.  READ  CONFIGURATION  FILE 

(pc  a  b  tc  cdll  cdl2  h31  es33  tan8  Qmech  psub  Esub  vsub  tsub  zO )  :=  READPRN(CC3L_config) 

a  =0.1016  b  =0.0625  tc  =0.00254  vc:  =  -^-^  vc  =0.485282 

cdll 

v  :=  0.5( vc -h  vsub)  v  =0.427641 

II.  MODE  SOLUTION 

In  this  section,  the  Rayleigh-Ritz  technique  is  used  to  find  the  fundamental  eigenvalue  and  eigenfunction  for 
an  8th  order  polynomial.  (N  is  preset  to  8  below.) 

UU  :  =  U(N,b,v)  TT  :=T(N,b,v) 

gvals  :=genvals(UU,TT) 

gvecs  :=genvecs(UU)TT) 

sgcomb  :  =  rsort(stack(gvalsT,  gvecs) ,  o) 

fundval  :  =  J sgcomb0  Q  fundval  =4.078196 


III.  FIND  THE  COEFFICIENTS  FOR  THE  DEFLECTION  FUNCTION 
fundcoef  :  =  submatrix(  sgcomb,  1  ,N,0,0) 

CC  :  =  scaledcoef(  fundcoef,  N,b) 


The  coefficients  are  normalized  so  that  the  maximum  deflection  is  one. 


IV.  PLOT  THE  NORMALIZED  STRESSES 


eta  :  =  b ,b  -t-  0.002..  1 


V.  CALCULATE  THE  SHAPE  INTEGRALS 


Ike  :=CCT  TT-CC 

Ike  =  0.225435 

la  :=  1  -  b2 

la  =0.996094 

Iv  :=  CCT- VV(N,b) 

Iv  =0.61 1179 

Iq  :=  CCT  QQ(N,b) 

Iq  =  1.076708 

Ipsum  :  =  CCT-Usum(N,b,v)-CC 

Ipsum  =4.412892 

Ipl  :=  CCT  Upl(N,b,  v)  CC 

Ipl  =1.1593 

IpO  :  =  Ipsum  -  Ipl 

IpO  =3.253593 

The  maximum  strain  is  at  the  inner  radius.  The  second  derivative  of  the  normalized  deflection  is 
calculated  below  (maxcurve)  and  the  peak  deflection  (i.e.,  deflection  at  r  =  a)  corresponding  to  a  strain  of 
one  is  saved  (wnmax). 

2 

"T 

maxcurve  :=CC  -SM(N,b)  wnmax  = -  wnmax  =  0.230222 

(tc+  zO)-maxcurve0 


Write  shape-integral  file: 

WRITEPRN(CC3L_defiect)  :=  (lp0Q  IplQ  Iq0  IvQ  IkeQ  la  wnmax  N  CC) 

Because  several  of  the  shape  integrals  are  calculated  through  matrix  operations,  they  remain  as 
matrices  even  though  having  dimension  1x1 .  The  subscript  index  is  used  to  convert  to  scalar  form  before 
writing.  The  deflection  coefficients  are  written  as  a  matrix  (Nxl)  as  CC. 


What  follows  are  predefinitions... 


N=8 


scaledcoef(  vect ,  N,  b )  = 


peakdefl*—  0 
for  iie  1..N 


incr*-  1  -  (ii-h  1  )*bu -h  ii*bll_r  1 
peakdefl*-  peakdefl  -t-  incrvect.._  J 


vectscaled< 

vectscaled 


vect 

peakdefl 


Usum(N,b,v)  s 


for  me  1 ..  N 
for  n  e  1 ..  N 

l_bm+n 

ul«-(m-h  1  )-(n  +■  1) - 

mtn 


u2*-  -  (n  +  1  ).bm-—  H-  -  (m  -  1  )-bn-^- 
n  m 

uam  _  i  n  _  t  (m+-  l)  (n-t-  1 )- (ul  +  u2+ -bm't'n-ln(b)) 


ua 


Upl(N,b,v)s 


for  m€  1..N 
for  ne  1..N 


ubb 

m 


1  ,n- 


n-(n-h  !)• 


l-bm^n 

mtn 


ub<-  ubb  -i-  ubbT 
ub 


U(N,b,v)  = 


uu<-  Usum(N,b,v)-  (1  -  v)-Upl(N,b,v) 


uu 


SM(N,b)  = 


for  me  1  ..N 
sm  l)*bm_  1 

m-  1  v  7 


sm 


T(N,b,v)  = 


for  me  1  ..N 
for  ne  1..N 


tU 


1  -  b1 


m  -4114-4 


m+-  n  -  4 


--  (n+  l)-b 


n  1  "  b 


m  4-4 


m+  4 


m-t-  3 


,  ,  n-t-l  1  -  b,,ITJ  ,  ...J-b1 

t2<-n-b  r - (m+-  l)b 


n-f-4 


n  +  4 


t3<—  (m  +■  1  )-(n  +  1  )  bm  — — -  (2'm-m-nn-n)-bm+n'1'1’- 


i  a  ■  m+1  1  — b"^3  ,m-(-n  +  2  1  — b 

t4<-m-b  r - -t-  rn-nb  T  T - 


m-  1  ,n_  1 


n  4-  3 
- 1 1 4- 12  4- 13  -T- 14 


tt 


VV(N,b)  = 


for  me  1  ..N 
matv 


m  -  1 


l-bm  +  3  ,  n,ml-b3 

- (m-  l)  b - 

m  4-  3  3 


4-  m-b1 


m4- 1  1  -  b 


2 -matv 


QQ(N,b)  = 


for  me  1  ..N 

matvm_  ^(m-h  l)*(l  -  bm) 


matv 


slope(  eta,  vect,  N,b)  = 


slopes- 0 
for  iie  1 ..  N 

slope*-  slope  4-  vect.._  j-(ii  H-  1  )-(etan-  b11) 
slope 


curve(  eta  ,  vect ,  N ,  b )  s 


curve*—  0 
for  iie  1  ..N 

curve*—  curve  4-  vect.._  j*(ii  +-  1  )-ii  eta' 
curve 


ii-  1 


Radial  Stress(  eta  5  vect ,  N ,  b ,  v)  s  curve(  eta ,  vect ,  N,  b )  4- 
TangentialStress( eta , vect ,N,b,v)=v- curve( eta,  vect , N, b )  4- 


v-slope(eta,  vect,N,b) 


eta 

slope(  eta,  vect,  N,b) 
eta 


PfN\SfATt; 
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Module  3:  Equivalent-Circuit  Parameters  --  Trilaminar  Flex  Disks 

This  worksheet  reads  the  configuration  and  shape-integral  files  for  a  particular  case  and  computes  the 
equivalent-circuit  parameters. 


V 

avg 


Alternate  form: 


In  these  equivalent  circuits,  the  mechanical  quantities,  force  (pA  or  acoustic  pressure  times  effective  area) 
and  velocity,  v,  are  related  to  the  electrical  quantities,  voltage,  e,  and  current,  i.  CEB  is  the  blocked 
electrical  capacitance  and  RD  is  the  dielectric  loss  (tan8/co*CEB).  The  mechanical  mass,  stiffness,  and 
mechanical  damping  are  m,  k,  and  Rm  respectively.  The  mechanical  resistance  is  computed  based  on  the 
mechanical  Q  of  the  ceramic  and  the  fraction  of  strain  energy  that  is  stored  in  the  ceramic  layer  of  the 
composite.  The  radiation  impedance,  Zrad,  is  represented  as  a  parallel  combination  of  resistance  and 
mass  (see  below).  Finally,  the  electromechanical  transduction  factor  (the  "turns  ratio")  is  <(>. 


I.  READ  FILES 
Read  configuration  file: 

(pc  a  b  tc  edit  cdl2  h31  es33  tanS  Qmech  psub  Esub  vsub  tsub  zO )  :=  READPRN(CC3L_config) 
a  =0.102  b  =0.063  tc  =  2.54- 10~3 


Read  shape-integral  file: 


(IpO  Ipl  Iq  Iv  Ike  la  wnmax  N  CC)  :  =  READPRN(CC3L_deflect) 

i  i  ~  zO  zO  Iz0\2 

acl  :=  1  -t-  2 —  ac2.-l  +  3 —  +  3*  — 

tc  tc  \  tc  / 

Aeff  :  =  7E*a2  la  Aeff  =  0.032 

II.  CALCULATE  EQUIVALENT-CIRCUIT  PARAMETERS 
Equivalent  mass: 
i  la2 

mass  :  =  4*7i*a  -Ike - (pc-tc-h  psub-tsub)  mass  =2.338 

Iv2 


zO 

as2  :=  1  -  3*-—-  -t-  3- 
tsub 


zO  \ 
tsub/ 


Equivalent  stiffness: 


o  ac  1 

fll  :=  cdl  1  -  0.75*es33*h31  •- — 

ac2 


fl2  :=cdl2-  0.75-es33-h312-— 

ac2 


csl  1 


Esub 
1  -  vsub2 


cs!2  :=  vsub-csll 


Iff  :  =  fl  1  ■  IpO  +■  f  1 2-  Ip  1 


Iss  =  csl  1  IpO h-  csl 2* Ipl 


stiffness  :  = 


3  a  Iv 


- (tc3-ac2-Iff  +■  tsub3*as2*Iss) 


stiffness  =2.056*  107 


Transduction  factor  for  parallel  connection  of  ceramic  layers: 


-  2- ( 7t-h3 1  es33-Iq  — tc-acl 

Iv 


(j)  =0.828 


Blocked  electrical  capacitance  for  parallel  connection  of  ceramic  layers: 
2 


CF.B  :=  2  g-s---'^1  a  -la 
tc 


CEB  =2.01*10 


-7 


Coupling  factor  (for  parallel  connection  of  ceramic  layers): 


kfl  :=  stiffness 


hr 

la2 


2  3  Iq2  acl2 
kf2  :=2-7l*8s33*h312Tc 

la  2 


k2  :  = 


kf2 


k2  =0.142 


kfl  -h  kf2 


Calculate  equivalent  Q  for  composite  structure: 

Qeff  :-Qmech-(  1  -t-  Uratio) 


Uratio  : 


tsub\3  as2  Iss 
to  /  ac2  Iff 


Rmech  :  = 


stiffness  mass 
4  mass  Qeff 


Calculate  circuit  elements  for  parallel  form  of  radiation  impedance.  This  form  is  preferred  over  the 
series  form  since,  in  the  parallel  form,  the  resistance  element  is  frequency  independent.  Furthermore,  the 
parallel  form  degrades  more  gracefully  for  ka  approaching  (or  greater  than)  one.  Also,  the  assumption  is 
made  that  these  transducer  elements  would  always  be  used  in  pairs  placed  back-to-back  and  driven  in 
phase.  In  this  drive  mode,  the  radiation  impedance  is  equivalent  to  that  of  a  single  element  in  an  infinite 
rigid  baffle. 


p_water  :=  1000  c_water  :  =  1500 


8  (  2\ 
m_rad  :  = - p_water-a- (jt-a  J*Ia 

3*71 


R_rad  :  =  1 .44  •  pjwater*  c_water 


Write  equivalent-circuit  parameters  to  a  file: 

WRlTEPRN(CC3L_equivckt)  :=(mass  stiffness  $  CEB  k2  tan5  Rmech  Rjrad  m_rad  Aeff) 


III.  MISCELLANEOUS  CALCULATIONS 

Unloaded  resonance  frequency: 

.  1  stiffness  , 

mech_res  :  = - -  mech_res  =471.953 

2-7C  4  mass 


Rough  estimate  of  water  loading  (piston  with  uniform  face  velocity): 


8  .  3  T 

rad_mass :  = - p_water* 7t- a  -la 

3*71 


co2water  :  = 


stiffness 

mass  +-  rad_mass 


pjwatera 

pc-2-tc 


=  2.632 


mech_res_water  :  = 


-W 


tt)2water 


mech_res_water  =  3 1 8.803 


2'K-  mech_res_water 
wa venum  ' = - 


c_water 


ka  : r  wavenum- a  ka  =0.136 


2  ks 

radR:-p  water-c  water-ft-a  -la -  radR  =445.971 

2 

Estimated  maximum  source  level  (if  strain  limited): 

[Note:  Maximum  strain  is  not  an  appropriate  failure  criterion  for  a  ceramic  material.  It  is  used  in  this 
draft  version  as  a  crude  estimate  of  limiting  performance.  Eventually,  the  more  appropriate  maximum 
tensile  stress  criterion  will  be  implemented.  The  often-used  von  Mises  stress  criterion  is  also  inappropriate 
for  failure  in  ceramic.] 


Enter  maximum  allowable  strain  Smax  :=  0.001 


Qmax  := 71* a2  - 2 -  71* mech  jres_water- wnmax- Smax- 1 v 

pi  :=  p_water-mech_res_water*Qrnax  SL  :=20-log(pl)  -}-  120 

Maximum  pressure  at  1  meter  (Pa):  pi  =  2.914*  103 

Maximum  SL  in  dB  with  respect  to  1  micropascal  at  one  meter:  SL  =  189.3 

Estimated  peak  TVR: 

What  follows  is  an  optimistic  estimate  of  the  peak  transmitting  voltage  response;  mechanical  loss 
in  the  ceramic  is  ignored. 

2  T 

tvr_peak  :=  p__water-mech_res_water-7i  a  Ta - 

radR 


TVR  in  pascals  at  one  meter  per  volt  input:  tvr_peak  =  19.122 

TVR  in  dB  with  respect  to  1  micropascal  at  one  meter  per  volt:  201og(  |tvr_peak| )  +  120  =  145.6 


Low-frequency  receive  response: 
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Module  1.ES2L  confia:  Configuration  Definition  -  Edge-Supported  Disk 

This  worksheet  is  specific  to  the  bilaminar  edge-supported  disk.  First,  the  physical  dimensions  are 
specified  and  the  piezoelectric  properties  file  is  read.  Then,  the  two-dimensional  properties  are  derived. 
Finally,  the  results  are  written  to  a  configuration  file. 


HINGED  EDGE  •  ts 

SUPPORT 


/.  SELECT  MATERIALS  [Highlighted  values  are  user  inputs] 

Choose  from  PZT4,  PZT5H,  or  PZT8  for  file  name  in  pzt  READPRN. 
pzt :  =  READPRN(  PZT4 ) 


sell  :=P*0f0-10' 

d31  :=pzt  lO-12— — 

°’3  volt 

Choose  from  steel,  aluminum,  or 
substrate  :=  READPRN( brass) 


se!2  :=pztQ  j*  10 


12  J_ 
Pa 


et33  =pztQ  2-s0 


P:  =  PZt0,4^ 

m 


tanS  :  =  pzt 


0,5 


brass  for  file  name  in  substrate  READPRN. 


Qmech  :=pztQ  6 


psub  :  =  substrateQ  Q — ~ 
m3 

psub  =  8.5*  103  *kg*m 


Esub  :=substrateQ  1-109-Pa  vsub  :-substrateQ  2 

Esub  =  1.04*  10U  *kg*m  1  *sec 


vsub  =0.37 


II.  SELECT  DIMENSIONS 


Enter  dimensions  of  ceramic  disk  (a  =  outside  radius;  tc  =  thickness) 
and  thickness  (tsub)  of  substrate  disk: 

[Enter  units  for  in,  mm,  cm,  or  m;  if  no  units  are  entered,  meters  will  be  assumed] 
a  :=  10cm  tc:=0.1-mm  tsub:=l-cm  b  :-0 


a  =  0.1*m  b  =  0 


tc  =  1*  10  4  •m 


III .  DERIVE  PROPERTIES  FOR  TWO-DIMENSIONAL  ANALYSIS 


sdll  : =  se 1 1  - 


d31^ 

et33 


,,  d312 

sdl 2  sel2 - 

st33 


sdl2 

sdll 


0.485 


cdl  1 


sdll 

sdll2-  sdl22 


cdl2  := 


-  sdl 2 

sdll2  -  sdl22 


cdl2 
cdl  1 


=  0.485 


es33  :  =  et33  - 


2'd312 
sel  1  +-  sel2 


h31 


d31 

et33  (sdl  1  -hsdl2) 


IV.  CALCULATE  LOCATION  OF  NEUTRAL  PLANE 

Origin  of  z-coordinate  at  neutral  plane;  zO  gives  distance  from  neutral  plane  to  ceramic  substrate 
interface: 


fs  :  =  - 


Esub 
1  -  vsub2 


zO  = 


tsub2-fs-  tc2-cd  1 1 
2*(tsub-fs+-  tc*cdl  1) 


zO  =4.951*  10  3  *m 


WARNING:  If  zO  is  negative,  then  the  neutral  plane  is  in  ceramic  and  some  part 
of  the  ceramic  will  be  in  tension  under  hydrostatic  load  unless  prestressed. 

V.  WRITE  CONFIGURATION  FILE 

PRN  files  cannot  handle  dimensions  so  ail  quantities  are  written  (in  SI)  without  units. 


out 


0,4  ' 


.  cdl  1 
cO 


out 


0,1 


out, 


0,5 


a 

mO 

cdl2 

cO 


OUt0,2  :  =  b 


out, 


0,6 


h31 

hO 


out, 


0,3  ' 


tc 

mO 


out 


0,7 


ss33 

eO 


outQ  8  :=tan8  outQ  g  :  =  Qmech 


out, 


0,10  ' 


psub 

po 


out 


0,11 


.  Esub 
cO 


zO 

mO 


out 


0,12 


:=  vsub 


out, 


0,13 


tsub 

mO 


out 


0,14  ' 


WRITEPRN(  ES2L_config)  : = out 

The  configuration  file  contains  the  following  quantities  in  this  order: 

density  of  ceramic,  p 
outer  radius  of  disk,  a 

ratio  of  inner  radius  to  outer  radius,  b  (always  zero  here) 

thickness  of  ceramic,  tc 

cdll  for  ceramic 

cd12  for  ceramic 

h31  for  ceramic 

ss33  for  ceramic 

loss  tangent  for  ceramic,  tan5 

mechanical  Q  of  ceramic 

density  of  substrate,  psub 

modulus  of  substrate,  Esub 

poisson's  ratio  of  substrate,  vsub 

thickness  of  substrate,  tsub 

z  location  of  neutral  plane,  zO,  with  respect  to  bottom  of  ceramic 
Unit  normalization  quantities  (predefined): 

ke  ,  „  ,  volt  A_,  farad 

p0  =  i— m0  =  l-m  cO=l-Pa  h0=l -  eO-1 

3  mm 

m 


£0  =  8.85410 


.  12  farad 


m 
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Module  2.ES2L  deflect:  Edge-Supported  Disk 

This  worksheet  reads  the  configuration  file,  calculates  the  first  axisymmetric  mode  deflection  function, 
and  then  calculates  all  of  the  shape  integrals.  The  shape  integrals  are  written  to  a  file  for  subsequent  use. 
In  this  version,  the  polynomial  Rayleigh-Ritz  method  is  used  to  compute  the  mode  deflection  function.  An 
eighth-order  polynomial  is  used.  This  module  is  expecting  the  configuration  file  from  version  2  --  that  is,  the 
configuration  with  arbitrary  location  of  the  neutral  plane  and  specific  properties  for  the  substrate. 


/.  READ  CONFIGURATION  FILE 

(pc  a  b  tc  edit  cdl2  h31  es33  tan5  Qmech  psub  Esub  vsub  tsub  zO)  :-READPRN(ES2L_config) 

a  =0  1  b  =0  tc  =0.0001  vc:=-^^  vc  =0.485282 

cdl  1 

v  :  =  0.5-(vc+ vsub)  v=  0.427641 

II.  MODE  SOLUTION 

In  this  section,  the  Rayleigh-Ritz  technique  is  used  to  find  the  fundamental  eigenvalue  and  eigenfunction  for 
an  8th  order  polynomial.  (N  is  preset  to  8  below.) 

UU  :=U(N,b,v)  TT  :=T(N,b,v) 

gvals  :=  genvals(UU,TT) 

gvecs  :  =  genvecs(  UU ,  TT ) 

sgcomb  :  =  rsort  (stack  (gvalsT,  gvecs) ,  o) 

fundval  ^^sgcomt^  Q  fundval  =5.116155 


III.  FIND  THE  COEFFICIENTS  FOR  THE  DEFLECTION  FUNCTION 
fundcoef :  =  submatrix(  sgcomb,  1 ,  N,  0 , 0 ) 

CC  :  =  scaledcoef( fundcoef, N,b) 

The  coefficients  are  normalized  so  that  the  maximum  deflection  is  one. 


IV.  PLOT  THE  NORMALIZED  STRESSES 


eta  :  =  b,b  -h  0.002..  1 


V.  CALCULATE  THE  SHAPE  INTEGRALS 


Ike  :=CCT  TT  CC 

Ike  =  0.140926 

la  :=  1  -  b2 

la  =  1 

Iv  :=  CCT  VV(N,b) 

Iv  —  0.445882 

Iq  :=  CCT  QQ(N,b) 

Iq  =  1.401826 

Ipsum  :=  CCT  Usum(N,b, v)-CC 

Ipsum  =4.813486 

Ipl  :  =  CCT-Upl(N,b,v)-CC 

Ipl  =  1.965116 

IpO  :=  Ipsum  -  Ipl 

IpO  =2.84837 

The  maximum  strain  is  at  the  inner  radius.  The  second  derivative  of  the  normalized  deflection  is 
calculated  below  (maxcurve)  and  the  peak  deflection  (i.e.,  deflection  at  r  =  a)  corresponding  to  a  strain  of 
one  is  saved  (wnmax). 

2 

maxcurve  :  =  CCXSM(N,b)  wnmax  :  = - - -  wnmax  =  0.732203 

(tc-h  z0)-maxcurveQ 


Write  shape-integral  file: 

WRITEPRN(ES2L_deflect)  :=  (lp0Q  IplQ  Iq„  Iv„  Ike0  la  wnmax  N  CC) 

Because  several  of  the  shape  integrals  are  calculated  through  matrix  operations,  they  remain  as 
matrices  even  though  having  dimension  1x1 .  The  subscript  index  is  used  to  convert  to  scalar  form  before 
writing.  The  coefficients  are  written  as  a  matrix  (Nxl)  as  CC. 


What  follows  are  predefinitions... 


N=8 


scaledcoef(  vect ,  N,  b )  = 


peakdefl*- { 
for  iie  1 .. 
peakdefl  <- 


vectscaled- 


vectscaled 


Usum(N,b,v)  = 


for  me  1  ..N 
for  ne  1  ..N 


ua 


t,n-  1 


ua 


Upl(N,b,  v)  = 


for  me  1 ..  N 
for  n  e  1 ..  N 


ubb 


n 

m  -  I  ,n-  1 


ub<  ubb  t-  ubbT 
ub 


U(N,b,  v)  = 


uu<-  Usum(N,b,  v)  -  ( 
uu 


SM(N,b)s 


for  me  2..N 

sm  ,<-0 
m  -  1 


sm 


2 


N 

-peakdefl  vect.._  ] 
vect 

peakdefl 


m  +  1  )2-(n  j-  1  )2 
m  -  n 


•(n  +  l)-(m-h  1) 
min 


1  -  v)-Upl(N,b,v) 


sm 


T(N,b,  v)  = 


for  me  1  ..N 
for  ne  1  ..N 


tt 


m-  1  ,n-  1 


_J _ 

m-t-  n-j-  4 


1  1  1 

- - 1 — 

m+-  3  n  -i-  3  2 


tt 


VV(N,b)  = 


for  me  1  ..N 


matv  —  1  — 

m  -  1 


2 

m-h  3 


matv 


QQ(N,b)n 


for  me  1 .. N 
matv  (mi  1) 

m-  1  v 


matv 


slope(eta,vect,N,b)  = 


slopes-  0 
for  iie  1 ..  N 

slope<-  slope  i-  vect.._  ^(ii -t  l  )-etan 
slope 


curve(  eta ,  vect ,  N ,  b )  ^ 


curvet—  0 
for  iie  1 .. N 

curves- curve -h  vect.. _  ^(ii 1 )- ii-eta11  1 
curve 


RadialStress(  eta ,  vect ,  N ,  b ,  v)  =  curve(  eta ,  vect ,  N,  b )  -t- 


v-  slope(  eta ,  vect ,  N,  b ) 
eta 


TangentialStress(  eta,  vect  ,N,b,v)  =  v*  curve(  eta ,  vect ,  N,  b )  +- 


slope(  eta,  vect,  N,b) 
eta 
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Flexural-Disk  Transducer  Analytical  Model  -  DRAFT  VERSION  3.0 

Module  3:  Equivalent-Circuit  -  Bilaminar,  Edge-Supported  Flex  Disks 

This  worksheet  reads  the  configuration  and  shape-integral  files  for  a  particular  case  and  computes  the 
equivalent-circuit  parameters. 


V 

avg 


Alternate  form: 


In  these  equivalent  circuits,  the  mechanical  quantities,  force  (pA  or  acoustic  pressure  times  effective  area) 
and  velocity,  v,  are  related  to  the  electrical  quantities,  voltage,  e,  and  current,  i.  CEB  is  the  blocked 
electrical  capacitance  and  RD  is  the  dielectric  loss  (tan§/oo*CEB).  The  mechanical  mass,  stiffness,  and 
mechanical  damping  are  m,  k,  and  Rm  respectively.  The  mechanical  resistance  is  computed  based  on  the 
mechanical  Q  of  the  ceramic  and  the  fraction  of  strain  energy  that  is  stored  in  the  ceramic  layer  of  the 
composite.  The  radiation  impedance,  Zrad,  is  represented  as  a  parallel  combination  of  resistance  and 
mass  (see  below).  Finally,  the  electromechanical  transduction  factor  (the  "turns  ratio")  is  <j>. 

I.  READ  FILES 
Read  configuration  file: 

(pc  a  b  tc  cdll  cdl2  h31  6s33  tan5  Qmech  psub  Esub  vsub  tsub  zO )  :  =  READPRN(ES2L_config) 

tc  =  1  *  1 0~4 


a  =0.1 


b  =0 


Read  shape-integral  file: 


(IpO  Ipl  Iq  Iv  Ike  la  wnmax  N  CC)  :-READPRN(ES2L_deflect) 


acl  :=  1  -h  2-—  ac2:=l  +  3-— -3-f— )  as2  :=  1  -  3  —  +  3-f— ) 

tc  tc  V  tc  /  tsub  \tsub/ 


Aeff  :=7t-a  -la 


Aeff  -  0.031 


//.  CALCULATE  EQUIVALENT-CIRCUIT  PARAMETERS 


Equivalent  mass: 


2  Iq*- 

mass  :  =  2*7t*a  -Ike - (pc-tc-t-  psub-tsub) 

Iv2 


mass  =  3.82 


Equivalent  stiffness: 


fll  :  =  cd  1 1  -  0.75-£s33-h31 


f  1 2  :  =  cd!2-  0.75-es33-h3r 


1  -  vsub 


cs!2  := vsub-csl 1 


Iff  :  =  fl  MpO-h  f  1 2-  Ip  1 


Iss  .=  csl  HpO  +■  cs  1 2* Ip  1 


stiffness  := AZL.ifL.  (tc3-ac2*Iff -h  tsub3-as2-Iss) 
~  2  T  2 
3  a  Iv 


stiffness  =  1. 16*  10 


Transduction  factor: 


7t-h31*Ss33Tq — tc-acl 
Iv 


0  =  1.461 


Blocked  electrical  capacitance: 


CEB  :=  es3^  ?■— la 


CEB  =2.482*10 


Coupling  factor: 


kfl  :=  stiffness- 


kf2  :  =  7i  £s33-h312-tc3  ^3_  — 
la  a2 


k2  := 


kf2 


kfl  +■  kf2 


k2  =7.352*  10 


Calculate  equivalent  Q  for  composite  structure: 


TT  .  I  tsub\3  as2  Iss 

Uratio  :=  [ - - 

tc  /  ac2  Iff 


Qeff  :=  Qmech-(  1  -t-  Uratio) 


„  ,  stiffness  mass 

Rmech  :=  I- - - — 

^  mass  Qeff 

Calculate  circuit  elements  for  parallel  form  of  radiation  impedance.  This  form  is  preferred  over  the 
series  form  since,  in  the  parallel  form,  the  resistance  element  is  frequency  independent.  Furthermore,  the 
parallel  form  degrades  more  gracefully  for  ka  approaching  (or  greater  than)  one.  Also,  the  assumption  is 
made  that  these  transducer  elements  would  always  be  used  in  pairs  placed  back-to-back  and  driven  in 
phase.  In  this  drive  mode,  the  radiation  impedance  is  equivalent  to  that  of  a  single  element  in  an  infinite 
rigid  baffle. 


1 

Tc 

I 


p_water  :=  1000  c_water  :=  1500 

m_rad  :  =  — ■ •p_water-a*(7i*a2)-Ia  R_rad  :  =  1.44-p_water-c_water*(7t*a2)*Ia 

3*71 

Write  equivalent-circuit  parameters  to  a  file: 

WRITEPRN(ES2L„equivckt)  :=(mass  stiffness  <j>  CEB  k2  tan8  Rmech  R_rad  m_rad  Aeff) 


ill .  MISCELLANEOUS  CALCULATIONS 

Unloaded  resonance  frequency: 

,  1  stiffness  ,  onn 

mech_res  :  = - -  mech_res  =  877.156 

2*71  ^ 


mass 


Rough  estimate  of  water  loading  (piston  with  uniform  face  velocity): 


rad_mass 


8  *  3  T 

= - p_water-7t-a  ‘la 

3-7t 


oo2  water :  = 


stiffness 

mass  +-  rad_mass 


p_water  a 


=  65.789 


pc-2-tc 


mech_res_water  := 


1 


4 


0)2water 


wavenum . 


2*7C 

2 -7i*mech_resjwater 
c  water 


mech_res_water  =  673.1 13 

ka  :  =  wavenum* a  ka  =0.282 


2  ka  3 

radR  :=p  water-c  water- 7t- a  -la -  radR  =  1.873*  10 

2 

Estimated  maximum  source  level  (if  strain  limited): 

[Note:  Maximum  strain  is  not  an  appropriate  failure  criterion  for  a  ceramic  material.  It  is  used  in  this 
draft  version  as  a  crude  estimate  of  limiting  performance.  Eventually,  the  more  appropriate  maximum 
tensile  stress  criterion  will  be  implemented.  The  often-used  von  Mises  stress  criterion  is  also  inappropriate 
for  failure  in  ceramic.] 


Enter  maximum  allowable  strain  Smax  0.001 


Qmax  : "  7i-  a2 •  2  -7i-  mech_res_water- wnmax-  Smax- 1 v 

pi  :=  |p_watermech_res_water-Qmax|  SL  :=20-log(pl)  +-  120 

Maximum  pressure  at  1  meter  (Pa):  pi  —  2.92*  104 

Maximum  SL  in  dB  with  respect  to  1  micropasca!  at  one  meter:  SL  =  209.3 

Estimated  peak  TVR: 

What  follows  is  an  optimistic  estimate  of  the  peak  transmitting  voltage  response;  mechanical  loss 
in  the  ceramic  is  ignored. 

2  T  $ 

tvr_peak  :=  p_water-mech_res_water-7t*a  Ta - 

radR 

TVR  in  pascals  at  one  meter  per  volt  input:  tvr„peak  =  16.489 

TVR  in  dB  with  respect  to  1  micropascal  at  one  meter  per  volt:  201og(  |tvr_peak| )  +  120  =  144.3 


Low-frequency  receive  response: 


T he  receiving  response  is  not  normally  important  but  can  be  useful  in  reciprocity  calculations. 


ffvs  := 


7ia2*Ia 


,  CEB*  stiffness 
1  + - 

<t>2 


FFVS  in  open-circuit  volts  per  pascal:  ffvs  =  1.581- 10 

FFVS  in  dB  with  respect  to  one  volt  per  micropascal:  201og(  |ffvs| )  -  120  =-196.0 

The  quantity  in  square  brackets  in  the  equation  for  ffvs  should  equal  02: 


- - - =7.352-10 

,  CEB- stiffness 

1  H - 

4> 


k2  =7.352-10 
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Module  1.ES3L  confia:  Configuration  Definition  -  Edge-Supported  Trilaminar  Disl 

This  worksheet  is  specific  to  the  three-layer  edge-supported  disk.  First,  the  physical  dimensions  are 
specified  and  the  piezoelectric  properties  file  is  read.  Then,  the  two-dimensional  properties  are  derived. 
Finally,  the  results  are  written  to  a  configuration  file. 


I 

mm 

HINGED  EDGE 
SUPPORT 


/.  SELECT  MATERIALS  [Highlighted  values  are  user  inputs] 

Choose  from  PZT4,  PZT5H,  or  PZT8  for  file  name  in  pzt  READPRN. 
pzt :  =  RE ADPRN(  PZT4 ) 


sell  :=pzt<>  o-10  •- 

d3l  :-pzt  ,-10"l2>- — 

°>J  volt 

Choose  from  steel,  aluminum, 
substrate  :  =  RE ADPRN( brass) 

k  e 

psub  :=substrateQ  Q-— 
rn 

psub  =  8.5*  103  -kg*m  3 


se!2  “pztQ  j-10'12—  st33  =pztQ  2*s0 


Pa 


p  =Pzto,ri 

m 


tanS  :  =  pzt 


0,5 


or  brass  for  file  name  in  substrate  READPRN. 


Qmech:=pzt0  6 


Esub  :=substrateQ  1*109-Pa  vsub  :  =  substrateQ  2 

Esub  =  1.04*1 011  ‘kg’m  l,sec 


vsub  =0.37 


//.  SELECT  DIMENSIONS 


Enter  dimensions  of  ceramic  disk  (a  =  outside  radius;  tc  =  thickness) 
and  thickness  (tsub)  of  substrate  disk: 

[Enter  units  for  in,  mm,  cm,  or  m;  if  no  units  are  entered,  meters  will  be  assumed] 


a  :=  10*  cm 
a  =  0.  I'm 


tc  :=  1  -  mm 


tsub  :=  1-cm 


b  =0 


tc  —  1  •  1 0  *m 


ts2 


tsub 


HI .  DERIVE  PROPERTIES  FOR  TWO-DIMENSIONAL  ANALYSIS 


sdll  :=sell 


d3T 

et33 


sd!2  :  =  se!2 - 


d3P 

et33 


sd  1 2 
sdll 


=  0.485 


cdl 1  :=- 


sdll 


sdll2-  sd!22 


cd!2  : 


-  sd  1 2 


sdll2-  sdi22 


cdl  2 
cdl  1 


=  0.485 


ss33  :=  et33  -  ■ 


2*d3 1 


sell  +-  se!2 


h31  :  =  - 


d3l 


et3 3 - ( sd  1 1  -h  sd!2) 


IV.  WRITE  CONFIGURATION  FILE 

PRN  files  cannot  handle  dimensions  so  all  quantities  are  written  (in  SI)  without  units. 


out, 


o,o ' 


out, 


0,4 


Q-l 

11 

outo,. 

a 

mO 

rxi 

O 

O 

b 

tc 

0UV>:=^ 

_  cdl  1 

outA  « 

.  _  cd  1 2 

°utn  , 

,_  h3 1 

es33 

°utn  7  "  . 

cO 


0,5 


cO 


hO 


eO 


outQ  g  :  =  tan5  outn  Q  =  Qmech  out. 


0,9 


0,10  * 


psub 

po 


out, 


0,11 


.  Esub 
cO 


OUt0, 12  vsub 


out, 


0,13 


isl 

mO 


out, 


0,14 


.  ts2 
mO 


WRITEPRN(  ES3L_config)  :=out 


The  configuration  file  contains  the  following  quantities  in  this  order: 


density  of  ceramic,  p 
outer  radius  of  disk,  a 

ratio  of  inner  radius  to  outer  radius,  b  (always  zero  here) 

thickness  of  ceramic,  tc 

cdll  for  ceramic 

cd12  for  ceramic 

h31  for  ceramic 

es33  for  ceramic 

loss  tangent  for  ceramic,  tanS 

mechanical  Q  of  ceramic 

density  of  substrate,  psub 

modulus  of  substrate,  Esub 

poisson's  ratio  of  substrate,  vsub 

half  thickness  of  substrate,  ts2 

location  of  neutral  plane  (same  as  half  substrate  thickness) 

Unit  normalization  quantities  (predefined): 

Vo  volt 

pO  =  l’—  mO=l*m  cO=l*Pa  hO=l -  e0  = 

3  m 


£0  =  8.854-10 


-12  farad 


farad 

m 


m 
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Module  2. ES3L  deflect:  Edge-Supported  Trilaminar  Disk 

This  worksheet  reads  the  configuration  file,  calculates  the  first  axisymmetric  mode  deflection  function, 
and  then  calculates  all  of  the  shape  integrals.  The  shape  integrals  are  written  to  a  file  for  subsequent  use. 
In  this  version,  the  polynomial  Rayleigh-Ritz  method  is  used  to  compute  the  mode  deflection  function.  An 
eighth-order  polynomial  is  used.  This  module  is  expecting  the  configuration  file  from  version  2  --  that  is,  the 
configuration  with  arbitrary  location  of  the  neutral  plane  and  specific  properties  for  the  substrate. 


/.  READ  CONFIGURATION  FILE 

(pc  a  b  tc  edit  cdl2  h31  ss33  tanS  Qmech  psub  Esub  vsub  tsub  zO)  :=READPRN(ES3L_config) 

a  =0.1  b=0  tc  =  0.001  vc:=-^^  vc  =0.485282 

cdl  1 

v  :=  0.5*(vc  +  vsub)  v=  0,427641 
II.  MODE  SOLUTION 

In  this  section,  the  Rayleigh-Ritz  technique  is  used  to  find  the  fundamental  eigenvalue  and  eigenfunction  for 
an  8th  order  polynomial.  (N  is  preset  to  8  below.) 

UU  :=  U(N,b,v)  TT  :=T(N,b,v) 

gvals  :=  genvals(UU,TT) 

g vecs  :  =  genvecs(  UU ,  TT ) 

sgcomb  :=  rsort (stack (gvalsT,  gvecs) ,  o) 

fundval  :  =  J sgcomb0  Q  fundval  =5.116155 


III.  FIND  THE  COEFFICIENTS  FOR  THE  DEFLECTION  FUNCTION 
fundcoef  :=  submatrix(  sgcomb,  1  ,N,0,0) 

CC  :  =  scaledcoef(  fundcoef,  N,b) 


The  coefficients  are  normalized  so  that  the  maximum  deflection  is  one. 


IV.  PLOT  THE  NORMALIZED  STRESSES 


eta  :=b,b  +-  0.002..  1 


V.  CALCULATE  THE  SHAPE  INTEGRALS 


Ike  :=CCT  TT-CC 

Ike  =  0.140926 

la  :=  1  -  b2 

la  =  1 

Iv  :=  CCT- VV(N,b) 

Iv=  0.445882 

Iq  :=  CCT  QQ(N,b) 

Iq  =  1.401826 

Ipsum  :  =  CCT-Usum(N,b,v)  CC 

Ipsum  =4.813486 

Ipl  =CCTUpl(N,b,  v)CC 

Ipl  =  1.965116 

IpO  :=  Ipsum  -  Ipl 

IpO  =2.84837 

The  maximum  strain  is  at  the  inner  radius.  The  second  derivative  of  the  normalized  deflection  is 
calculated  below  (maxcurve)  and  the  peak  deflection  (i.e.,  deflection  at  r  =  a)  corresponding  to  a  strain  of 
one  is  saved  (wnmax). 

2 

maxcurve  :=  CCT-SM(N,b)  wnmax  :  = -  wnmax  =0.616393 

(tc-h  zO)*maxcurve0 


Write  shape-integral  file: 

WRITEPRN(ES3L_deflect)  -(lp00  IplQ  IqQ  IvQ  IkeQ  la  wnmax  N  CC) 

Because  several  of  the  shape  integrals  are  calculated  through  matrix  operations,  they  remain  as 
matrices  even  though  having  dimension  1x1 .  The  subscript  index  is  used  to  convert  to  scalar  form  before 
writing.  The  coefficients  are  written  as  a  matrix  (Nxl)  as  CC. 


What  follows  are  predefinitions... 


N=8 


scaledcoef(  vect ,  N ,  b )  = 


peakdefl*- 0 
for  lie  1..N 

peakdefl*-  peakdefl  -t-  vect.._  { 


,  ,  vect 

vectscaled* - 

peakdefl 

vectscaled 


Usum(N,b,  v)  = 


for  me  1  ..N 
for  ne  1  ..N 


ua 


m  -  !  ,n -  1 


(m  +-  1  )2-(n  j-  1  )2 
m  -  n 


ua 


Upl(N,b,  v)  = 


for  me  1  ..N 
for  ne  1 .. N 


ubb  .  .* 

m  -  1 ,  n  -  1 
T 

ub<-  ubb  i-  ubb 
ub 


n-(n-j-  1)‘(ith-  1) 
m  -i  -  n 


U(N,b,  v)  = 


uu<—  Usum(N,b,  v)  -  ( l  -  v)  Upl(N,b,  v) 
uu 


SM(N,b)  = 


for  me  2..N 

sm  ,<-0 
m  -  i 


sm 


2 


sm 


T(N,b,v)s 


for  me  1 .. N 
for  n  e  1 ..  N 


m  l,n  1  mtn-h4 


1  1  1 

- +  _ 

m-h  3  nt  3  2 


tt 


VV(N,b)  = 


for  me  1  ,.N 


matv 

m  -  l 


1  - 


2 

m-h  3 


matv 


QQ(N,b)s 


for  me  1 .. N 

matv  .(-fnn-l) 
m-  1  v  7 

matv 


slope( eta,  vect, N,b)  = 


slopes- 0 
for  iie  1 ..  N 


slopes-  slope  i-  vect.._  j  (ii  i-  1  )’eta11 
slope 


curve  ( eta ,  vect ,  N ,  b )  = 


curve  0 
for  iie  1  ..N 

curve*— curve -h  vect..  ,*(ii *h  1  )* ii* eta11  1 
n  —  1  v  ' 

curve 


RadialStress(  eta ,  vect ,  N ,  b ,  v)  =  curve(  eta ,  vect ,  N ,  b )  -h 


V'  slope(  eta ,  vect ,  N ,  b ) 
eta 


T  angentialStress(  eta ,  vect ,  N ,  b ,  v)  =  v*  curve(  eta ,  vect ,  N ,  b  )  -h 


slope(  eta,  vect,  N,b) 
eta 
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Module  3:  Equivalent-Circuit  -  Trilaminar,  Edge-Supported  Flex  Disks 

This  worksheet  reads  the  configuration  and  shape-integral  files  for  a  particular  case  and  computes  the 
equivalent-circuit  parameters. 


V 

avg 


Alternate  form: 


In  these  equivalent  circuits,  the  mechanical  quantities,  force  (pA  or  acoustic  pressure  times  effective  area) 
and  velocity,  v,  are  related  to  the  electrical  quantities,  voltage,  e,  and  current,  i.  CEB  is  the  blocked 
electrical  capacitance  and  RD  is  the  dielectric  loss  (tan8/af  CEB).  The  mechanical  mass,  stiffness,  and 
mechanical  damping  are  m,  k,  and  Rm  respectively.  The  mechanical  resistance  is  computed  based  on  the 
mechanical  Q  of  the  ceramic  and  the  fraction  of  strain  energy  that  is  stored  in  the  ceramic  layer  of  the 
composite.  The  radiation  impedance,  Zrad,  is  represented  as  a  parallel  combination  of  resistance  and 
mass  (see  below).  Finally,  the  electromechanical  transduction  factor  (the  "turns  ratio”)  is  <j>. 

I.  READ  FILES 
Read  configuration  file: 

(pc  a  b  tc  edit  cdl2  h31  es33  tanS  Qmech  psub  Esub  vsub  tsub  zO )  :=READPRN(ES3L_config) 
a  =0.1  b  =0  tc  =  1*1 0_ 3 


Read  shape-integral  file: 


(IpO  Ipl  Iq  Iv  Ike  la  wnmax  N  CC)  :  =  READPRN(ES3L_deflect) 

\  2  _ a  /  \  2 


acl  :=  1  +  2 


zO 

tc 


*2  :=  It  3— +  3.(5? 

tc  \tc 


*  1  <■>  zO  /  zO  \ 

as2  :=  1  -  3 - +-  3*  - 

tsub  Itsub/ 


Aeff  :  =  7i*a  -la 


Aeff  =  0.031 


II.  CALCULATE  EQUIVALENT-CIRCUIT  PARAMETERS 
Equivalent  mass: 

2  la2 

mass  :  =  4-7t*a  -Ike - (pc-tc-t-  psub-tsub)  mass  =4.463 

Iv2 

Equivalent  stiffness: 


fll  cdl  1  -  0.75-ss33-h3l 


2  acT 
ac2 


f!2  :=cdl2~  0.75-ss33-h31 


2  acl 
ac2 


csl  1  :  =  - 


Esub 
1  -  vsub2 


cs!2  :=  vsub-csl  1 


Iff  :=  fl  MpO -{-  fl2*Ipl 


Iss  =csl  1-IpO-f-  cs  1 2*  Ip  1 


stiffness  (tc3*ac2’Iff-i-  tsub3-as2Tss) 

2  T  2  ' 

3-a  Iv 


stiffness  =  1.747*  10 


Transduction  factor  for  parallel  connection  of  ceramic  layers: 


:  =  -2*  (xch31  ■es33Tq-— *tc-acl 
Iv 


=  3.213 


Blocked  electrical  capacitance: 


CEB  :=  2’(£s33'7I'a  )..Ia 
tc 


CEB  =4.965-10 


Coupling  factor  for  parallel  connection  of  ceramic  layers: 


Iv 

kfl  :  =  stiffness - 

la2 


kf2  :  =  2-7t*ss33*h312-tc3*— 

I*  a2 


k2  := 


kf2 


k2  =0.106 


kfl  -h  kf2 


Calculate  equivalent  Q  for  composite  structure: 

Uratio  :=  —  Qeff  :=Qmech  (l  +  Uratio) 

\  tc  /  ac2  Iff 


^  ,  stiffness  mass 

Rmecn  := - 

mass  Qeff 

Calculate  circuit  elements  for  parallel  form  of  radiation  impedance.  This  form  is  preferred  over  the 
series  form  since,  in  the  parallel  form,  the  resistance  element  is  frequency  independent.  Furthermore,  the 
parallel  form  degrades  more  gracefully  for  ka  approaching  (or  greater  than)  one.  Also,  the  assumption  is 
made  that  these  transducer  elements  would  always  be  used  in  pairs  placed  back-to-back  and  driven  in 
phase.  In  this  drive  mode,  the  radiation  impedance  is  equivalent  to  that  of  a  single  element  in  an  infinite 
rigid  baffle. 


p_water  :  =  1 000  cjwater :  =  1 500 

m_rad  'p_water-a-(7i-a2)-Ia  R_rad  :=  l.44*p_water*c_waten (n-a2)-Ia 

3*71 


Write  equivalent-circuit  parameters  to  a  file: 


WRITEPRN(  ES3L_equivckt)  :  =  ( mass  stiffness  cj)  CEB  k2  tan8  Rmech  R_rad  m_rad  Aeff) 


III.  MISCELLANEOUS  CALCULATIONS 


Unloaded  resonance  frequency: 


mech_res  :  = - 

2-Jt 


stiffness 


mech_res  =995.732 


mass 


Rough  estimate  of  water  loading  (piston  with  uniform  face  velocity): 


racLmass 


8  *  3  T 

;= - p_water-7t*a  -la 

3tc 


co2water :  = 


stiffness 

mass  +  rad_mass 


p_water-a 

pc-2*tc 


=  6.579 


mech_res_water  :  =  o2water 
2-n 


2  •  7i*  mech_res_water 

wavenum  := - 

c  water 


mech_res_water  =  787.8 


ka  :=  wavenum-a  ka=0.33 


2  J(2  3 

radR:  =  p  water-c  water-Tt-a  -la -  radR  =2.566*  10 

2 

Estimated  maximum  source  level  (if  strain  limited): 

[Note:  Maximum  strain  is  not  an  appropriate  failure  criterion  for  a  ceramic  material.  It  is  used  in  this 
draft  version  as  a  crude  estimate  of  limiting  performance.  Eventually,  the  more  appropriate  maximum 
tensile  stress  criterion  will  be  implemented.  The  often-used  von  Mises  stress  criterion  is  also  inappropriate 
for  failure  in  ceramic.] 


Enter  maximum  allowable  strain  Smax  :  =  0.001 


2 

Qmax  :~7t-a  •2*rc*mechjres_waterwnmax*SmaxTv 
p  1  :  =  |  p_water mech_res_water- Qmax |  SL  :  =  20* log(  pi )  +  120 

Maximum  pressure  at  1  meter  (Pa):  pi  =  3.367*  104 

Maximum  SL  in  dB  with  respect  to  1  micropascal  at  one  meter:  SL  =210.5 

Estimated  peak  TVR: 

What  follows  is  an  optimistic  estimate  of  the  peak  transmitting  voltage  response;  mechanical  loss 
in  the  ceramic  is  ignored. 


tvr_peak  :  =  p_water  mech_res_water-7i-  a  •  la* 


4> 

radR 


TVR  in  pascals  at  one  meter  per  volt  input:  tvr_peak  =  30.989 

TVR  in  dB  with  respect  to  1  micropascal  at  one  meter  per  volt:  20-log(  |tvr_peak| )  +■  120  =  149.8 


Low-frequency  receive  response: 


The  receiving  response  is  not  normally  important  but  can  be  useful  in  reciprocity  calculations. 


ffvs  :  = 


7C*a2Ta  f 


1 


1  + 


CEB- stiffness 


FFVS  in  open-circuit  volts  per  pascal:  ffvs  =  1.04*10 

FFVS  in  dB  with  respect  to  one  volt  per  micropascal:  20-log(  |ffvs| )  -  120  =-179.7 

The  quantity  in  square  brackets  in  the  equation  for  ffvs  should  equal  02: 


.  . -  =0.106  k2  =0.106 

,  CEB- stiffness 
1  H- - 


2 
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Flexural-Disk  Transducer  Analytical  Model  -  DRAFT  VERSION  3.0 


Performance  Module 

This  worksheet  can  be  used  with  any  of  the  flexural-disk  configurations.  Several  performance-related 
measures  are  calculated  from  the  basic  equivalent-circuit  representation.  These  include:  (1)  unloaded 
admittance,  (2)  water-loaded  admittance,  (3)  transmitting  voltage  response,  and  (4)  transmitting  current 
response.  For  (3)  and  (4)  the  radiation  impedance  is  taken  to  be  that  of  a  piston  in  an  infinite,  rigid  baffle. 
This  is  appropriate  either  for  a  truly  baffled  transducer  or  for  a  transducer  that  was  constructed  from  TWO 
flexural-disk  elements  mounted  back-to-back  as  would  normally  be  done  in  a  high-performance  system. 
The  two-element  system  is  the  image  equivalent  of  a  single  projector  in  an  infinite,  rigid  baffle.  Remember 
that  the  electrical  characteristics  (e.g.,  (1)  and  (2)  above)  are  for  a  single  element.  If  a  two-element  system 
is  being  modeled,  the  user  must  adjust  the  electrical  parameters  according  to  the  electrical  connection  of 
the  elements  (series  or  parallel). 


Required  user  input  is  the  name  of  the  file  containing  the  equivalent-circuit  parameters: 


(mass  stiffness  <j)  CEB  k2  tan8  Rmech  R_rad  m_rad  Aeff )  :=  READPRN(CC2L_equivckt) 


Restore  physical  SI  units: 
mass  :  =  mass*mass_unit 
<j)  :  =  (|)'(j>_unit 

Rmech  :  =  Rmech* Rmech_un it 
m_rad  m__rad*mass_unit 


stiffness  :  =  stiffness- stiffness_unit 
CEB  :=  CEB*capacitance_unit 

R_rad  :  =  R_rad*Rmech„unit 
Aeff  :=  Aeff*area_unit 


If  you  want  to  model  a  closed,  fluid-filled  cavity  behind  the  flexural  disk,  set  'cavity'  equal  to  one;  if  not, 
make  sure  that  'cavity'  is  set  to  zero  and  skip  entry  of  cavity/fluid  properties. 


cavity  :=  0 


Enter  optional  fluid-cavity  properties  (set  'cavity' to  zero  to  disable 

Enter  volume,  vc,  of  cavity;  density,  pf,  of  cavity  fluid;  and  sound  speed,  cf,  of  cavity  fluid: 
[For  a  back-to-back  two  element  transducer,  enter  only  half  the  total  cavity  volume.] 

vc:=0.0016m3  pf:=800-^  cf:=  1400  — 

3  sec 


[Make  sure  to  enter  units  for  quantities;  if  no  units  are  entered,  SI  will  be  assumed] 


Cavity  stiffness: 


kcpfcf2.^! 


vc 


stiffness  =4.157*  107  ‘kg’sec”2 
kc  =  2.229*  1 08  *  kg*  sec  2 


stiffness  :=  stiffness  +  kc- cavity 


The  unloaded  resonance  frequency  and  Q  are 


^  _  1  stiffness 

2-n  4  mass 


fO  =  1.027*  103  -sec  1 


Q.unload  :=2-7t-fO-^5L  Q_unload  =  1.146- 103 
Rmech 


The  electrical  equivalents  of  the  mechanical  mass,  compliance,  and  resistance  are 

Rmech 


mass  e  : 


mass 


Cm_e  ;  =  - 


0 


stiffness 


Rmech__e  :  =  - 


The  admittance  calculation  is  broken  into  three  parts:  Y1  is  the  dielectric  loss,  Y2  is  the  blocked 
capacitance,  and  Y3  is  the  electrical  equivalent  of  the  mechanical  section. 


Yl(w)  :  =  w-CEB'tanS  Y2(w)  :=j  -w-CEB  Y3(w)  :  =  - 


Rmech_e-hj  -w*mass_e-t- 


1 


j  •  w-  Cm_e 


Y_unloaded( f)  :=  Yl(2-n-f)  +■  Y2(2-7i-f)  +•  Y3(2-*-f) 


The  awkward  part  of  plotting  the  unloaded  admittance  is  that  the  unloaded  Q  is  high  and  the  resonance 
peak  is  sharp.  The  fmat  routine  produces  a  matrix  (fx)  of  unevenly  spaced  frequencies  so  that  the 
resonance  is  well  defined  in  the  plot  without  excessive  calculation  of  points  away  from  the  resonance.  For 
water-loaded  calculations,  this  is  unnecessary. 

fx  : -  fmat( fO , Q_unload , k2 )  nn  :=  0..  length(fx)  -  1 

ymat(fx)  :=  j  for  mme  0..  length(fx)  -  1 

Wmm^Y— unloaded(fxmm)  Y_U  :=ymat(fx) 

1  yy 


UNLOADED  ADMITTANCE 


Frequency  [Hz] 


Calculate  the  water-loaded  resonance  frequency  and  the  Q: 


~  1  stiffness 

fOw := - 


2-71  v  mass  -t  m_rad 


xx 


/2-7t-fDw*m_rad 
\  R_rad 


R  at  fO  rad  R  rad 


(  xx 

1 1  +  XX 


Q_water  :  =  2’Tt-fO- 


mass  +  m„rad 
Rmech  h-  R_at_fD_rad 


fOw  =  745.664'  sec  Q_water  =22.728 


Calculate  the  electrical  equivalents  of  the  radiation  impedance  terms  and  calculate  the  electrical  equivalent 
of  the  radiation  impedance  (using  the  parallel  mass/resistance  model): 


,  m  rad 
mrad_e :  = - 


Rrad_e  :  =  — — ■ 


„  ,  ,  x  j  ‘W-Rrad  e-mrad  e 

Zrad_e(w)  - 


Rrad_e  +■  j  •  w*  mrad_e 


2 


2 


Calculate  the  admittance  term  for  the  mechanical  side  including  the  radiation  load: 


Y4(w)  := - - - 

Rmech_e  +-  j  •  w-  mass_e  4 - +-  Zrad_e(  w) 

j  -wCm_e 

Calculate  and  plot  the  total  admittance.  An  even  increment  in  frequency  is  suitable  since  the  Q  should  be 
much  lower  with  the  water  load. 


Y_water(0  :=  Yl(2-n-f)  +  Y2(2rc-f)  +  Y4(2-rcf) 

ffl  =0.5-f0w  ff2  :=  1.5-fDw  dff:=ff2~  —  ffx  :=  ffl ,ffl  4-  dff..  ff2 

200 


WATER-LOADED  ADMITTANCE 


Calculation  of  Transmitting  Voltage  Response 


The  pressure  (in  Pa)  at  one  meter  per  volt  drive  is 

rho  water  :=1000  —  p_per_volt(f)  :  =  Aeff- rho_water- f- Y4^  — 

3  (h 

m  Y 

or,  in  dB  with  respect  to  one  micropascal  at  one  meter  per  volt, 

dB_TVR_ref  :=  10"6-— •  1  m 
volt 

dB  TVR(f)  :=20  logfP=E5!=^2!!tO'| 
dB  TVR  ref 


TRANSMITTING  VOLTAGE  RESPONSE 


Frequency  [Hz] 


Calculation  of  Transmitting  Current  Response 

The  pressure  (in  Pa)  at  one  meter  per  ampere  drive  is 

_p_per_volt(f) 

p_per_amp(  f)  •- - - 

Y_water(  f) 

or,  in  dB  with  respect  to  one  micropascal  at  one  meter  per  ampere, 


dB_TCR_ref  :=  10'6— lm 


dB_TCR(f)  :=  20-log 


amp 

!  p_per_amp(f)\ 
\  dB_TCR_ref } 


1200 


TRANSMITTING  CURRENT  RESPONSE 


Predefinition  of  routine  to  calculate  frequencies  for  unloaded  admittance  plot: 


fmat(fO,Q,k2)  = 


fit- 0.5 
f2<—  0.6 


fdip< -  ■■■■ 

-  k2 


dpm<—  fdip  -  1 
f3<-  1  ±  0.375-dpm 
f4«-  1  +■  0.625-dpm 
f5<-  1.1 -fdip 
f6<—  1.2- fdip 
2 

fom<—  1 - 

Q 

fop< —  1  -h  — 

Q 


fdipm<-  fdip  -  — 
fdipp<—  fdip  -h  — 

Q 

f5<—  if(f5<1.4, 1.4,  f5) 
f6*-  if(f6<1.5, 1.5,f6) 
f2  <-  if(  f2  <fom ,  fom ,  f2 ) 
f3<-  if(f3>fop,fop,f3) 
f4*-  if(  f4<fdipm,  fdipm,  f4) 
f5  ^  if(  f5  >fdipp ,  fdipp ,  f5 ) 
f2  -  fl 


dfa 


100 


cum_m<—  0 
for  mme  0.. 99 


ff  <-fl-hmm-dfa 

mm 


dfb< 


O-  f2 
50 


cum_m<-  cum„m  4-  1 00 

for  mme  0.. 49 

ff  mm-dfb 

cum_m  -h  mm 


dfc 


f4  -  f3 
25 


cum_m<—  cum_m  h-  50 
for  mme  0 ..  24 


ff 


cum_m  4-  mm 


dfd< 


f5  -  f4 
50 


f3  4-  mm- dfc 


cum_m<-  cum_m  4-25 
for  mme  0..49 


ff 


cum_m  4-  mm 


■  f4  4-  mm-dfd 


dfe 


f 6-  f5 
100 


cum_m<—  cum_m  +-  50 
for  mme  0.. 99 


ff 

cum_m  4-  mm 

ff-fO 


-f5  4-  mm- dfe 


Predefinitions  of  unit  factors 

mass_unitsl*kg 


stiffness_unit=  1 


newton 


m 


(j)_unit=  1 


newton 


^  ,  newton- sec 

Rmech_unit  =  1  • - - 


m 


area  units  l*m 


capacitance_unit  *  1  -  farad 


volt 


Flexural-Disk  Transducer  Model 


T.  B.  Gabrielson 


Appendix  I.  MathCad  Materials  Files 
PZT4.pm 

12.3  -4.05  1300  -122  7600  0.004  500 

PZT5H.pm 

16.5  -4.78  3400  -274  7500  0.02  70 
PZT8.pm 

11.5  -3.7  1000  -97  7500  0.004  1050 

brass.pm 

8500  104  0.37 

aiuminum.pm 

2700  70  0.33 

steel.pm 

7700  195  0.28 
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